{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "# Plot style\n",
    "sns.set()\n",
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (4, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".pquote {\n",
       "  text-align: left;\n",
       "  margin: 40px 0 40px auto;\n",
       "  width: 70%;\n",
       "  font-size: 1.5em;\n",
       "  font-style: italic;\n",
       "  display: block;\n",
       "  line-height: 1.3em;\n",
       "  color: #5a75a7;\n",
       "  font-weight: 600;\n",
       "  border-left: 5px solid rgba(90, 117, 167, .1);\n",
       "  padding-left: 6px;\n",
       "}\n",
       ".notes {\n",
       "  font-style: italic;\n",
       "  display: block;\n",
       "  margin: 40px 10%;\n",
       "}\n",
       "img + em {\n",
       "  text-align: center;\n",
       "  display: block;\n",
       "  color: gray;\n",
       "  font-size: 0.9em;\n",
       "  font-weight: 600;\n",
       "}\n",
       "</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<style>\n",
    ".pquote {\n",
    "  text-align: left;\n",
    "  margin: 40px 0 40px auto;\n",
    "  width: 70%;\n",
    "  font-size: 1.5em;\n",
    "  font-style: italic;\n",
    "  display: block;\n",
    "  line-height: 1.3em;\n",
    "  color: #5a75a7;\n",
    "  font-weight: 600;\n",
    "  border-left: 5px solid rgba(90, 117, 167, .1);\n",
    "  padding-left: 6px;\n",
    "}\n",
    ".notes {\n",
    "  font-style: italic;\n",
    "  display: block;\n",
    "  margin: 40px 10%;\n",
    "}\n",
    "img + em {\n",
    "  text-align: center;\n",
    "  display: block;\n",
    "  color: gray;\n",
    "  font-size: 0.9em;\n",
    "  font-weight: 600;\n",
    "}\n",
    "</style>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\newcommand\\bs[1]{\\boldsymbol{#1}}\n",
    "\\newcommand\\norm[1]{\\left\\lVert#1\\right\\rVert}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span class='notes'>\n",
    "    This content is part of a series following the chapter 2 on linear algebra from the [Deep Learning Book](http://www.deeplearningbook.org/) by Goodfellow, I., Bengio, Y., and Courville, A. (2016). It aims to provide intuitions/drawings/python code on mathematical theories and is constructed as my understanding of these concepts. You can check the syllabus in the [introduction post](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-Introduction/).\n",
    "</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "We saw that not all matrices have an inverse. It is unfortunate because the inverse is used to solve system of equations. In some cases, a system of equation has no solution, and thus the inverse doesn't exist. However it can be useful to find a value that is almost a solution (in term of minimizing the error). We will see for instance how we can find the best-fit line of a set of data points with the pseudoinverse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.9 The Moore-Penrose Pseudoinverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Moore-Penrose pseudoinverse is a direct application of the SVD (see [2.8](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-2.8-Singular-Value-Decomposition/)). But before all, we have to remind that systems of equations can be expressed under the matrix form.\n",
    "\n",
    "As we have seen in [2.3](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-2.3-Identity-and-Inverse-Matrices/), the inverse of a matrix $\\bs{A}$ can be used to solve the equation $\\bs{Ax}=\\bs{b}$:\n",
    "\n",
    "$$\n",
    "\\bs{A}^{-1}\\bs{Ax}=\\bs{A}^{-1}\\bs{b}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\bs{I}_n\\bs{x}=\\bs{A}^{-1}\\bs{b}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\bs{x}=\\bs{A}^{-1}\\bs{b}\n",
    "$$\n",
    "\n",
    "But in the case where the set of equations have 0 or many solutions the inverse cannot be found and the equation cannot be solved. The pseudoinverse is $\\bs{A}^+$ such as:\n",
    "\n",
    "$$\n",
    "\\bs{A}\\bs{A}^+\\approx\\bs{I_n}\n",
    "$$\n",
    "\n",
    "minimizing\n",
    "\n",
    "$$\n",
    "\\norm{\\bs{A}\\bs{A}^+-\\bs{I_n}}_2\n",
    "$$\n",
    "\n",
    "The following formula can be used to find the pseudoinverse:\n",
    "\n",
    "$$\n",
    "\\bs{A}^+= \\bs{VD}^+\\bs{U}^T\n",
    "$$\n",
    "\n",
    "with $\\bs{U}$, $\\bs{D}$ and $\\bs{V}$ respectively the left singular vectors, the singular values and the right singular vectors of $\\bs{A}$ (see the SVD in [2.8](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-2.8-Singular-Value-Decomposition/)). $\\bs{A}^+$ is the pseudoinverse of $\\bs{A}$ and $\\bs{D}^+$ the pseudoinverse of $\\bs{D}$. We saw that $\\bs{D}$ is a diagonal matrix and thus $\\bs{D}^+$ can be calculated by taking the reciprocal of the non zero values of $\\bs{D}$.\n",
    "\n",
    "This is a bit crude but we will see some examples to clarify all of this.\n",
    "\n",
    "### Example 1.\n",
    "\n",
    "Let's see how to implement that. We will create a non square matrix $\\bs{A}$, calculate its singular value decomposition and its pseudoinverse.\n",
    "\n",
    "$$\n",
    "\\bs{A}=\\begin{bmatrix}\n",
    "    7 & 2\\\\\\\\\n",
    "    3 & 4\\\\\\\\\n",
    "    5 & 3\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.16666667, -0.10606061,  0.03030303],\n",
       "       [-0.16666667,  0.28787879,  0.06060606]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[7, 2], [3, 4], [5, 3]])\n",
    "U, D, V = np.linalg.svd(A)\n",
    "\n",
    "D_plus = np.zeros((A.shape[0], A.shape[1])).T\n",
    "D_plus[:D.shape[0], :D.shape[0]] = np.linalg.inv(np.diag(D))\n",
    "\n",
    "A_plus = V.T.dot(D_plus).dot(U.T)\n",
    "A_plus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now check with the `pinv()` function from Numpy that the pseudoinverse is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.16666667, -0.10606061,  0.03030303],\n",
       "       [-0.16666667,  0.28787879,  0.06060606]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.pinv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks good! We can now check that it is really the near inverse of $\\bs{A}$. Since we know that\n",
    "\n",
    "$$\\bs{A}^{-1}\\bs{A}=\\bs{I_n}$$\n",
    "\n",
    "with\n",
    "\n",
    "$$\\bs{I_2}=\\begin{bmatrix}\n",
    "    1 & 0 \\\\\\\\\n",
    "    0 & 1\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   2.63677968e-16],\n",
       "       [  5.55111512e-17,   1.00000000e+00]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_plus.dot(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not bad! This is almost the identity matrix!\n",
    "\n",
    "A difference with the real inverse is that $\\bs{A}^+\\bs{A}\\approx\\bs{I}$ but $\\bs{A}\\bs{A}^+\\neq\\bs{I}$.\n",
    "\n",
    "Another way of computing the pseudoinverse is to use this formula:\n",
    "\n",
    "$$\n",
    "(\\bs{A}^T\\bs{A})^{-1}\\bs{A}^T\n",
    "$$\n",
    "\n",
    "The result is less acurate than the SVD method and Numpy `pinv()` uses the SVD ([cf Numpy doc](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.pinv.html)). Here is an example from the same matrix $\\bs{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.16666667, -0.10606061,  0.03030303],\n",
       "       [-0.16666667,  0.28787879,  0.06060606]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_plus_1 = np.linalg.inv(A.T.dot(A)).dot(A.T)\n",
    "A_plus_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case the result is the same as with the SVD way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the pseudoinverse to solve a overdetermined system of linear equations\n",
    "\n",
    "In general there is no solution to overdetermined systems (see [2.4](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-2.4-Linear-Dependence-and-Span/) ; [Overdetermined systems](https://en.wikipedia.org/wiki/Overdetermined_system)). In the following picture, there is no point at the intersection of the three lines corresponding to three equations:\n",
    "\n",
    "<img src=\"images/overdetermined-system-equations-python.png\" width=\"300\" alt=\"Example of three linear equations in 2 dimensions: this is an overdetermined system\" title=\"Overdetermined system of equations\">\n",
    "<em>There is more equations (3) than unknowns (2) so this is an overdetermined system of equations</em>\n",
    "\n",
    "The pseudoinverse solve the system in the least square error perspective: it finds the solution that minimize the error. We will see this more explicitly with an example.\n",
    "\n",
    "<span class='pquote'>\n",
    "    The pseudoinverse solve the system in the least square error perspective\n",
    "</span>\n",
    "\n",
    "### Example 2.\n",
    "\n",
    "For this example we will consider this set of three equations with two unknowns:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "-2x_1 + 2 = x_2 \\\\\\\\\n",
    "4x_1 + 8 = x_2 \\\\\\\\\n",
    "-1x_1 + 2 = x_2\n",
    "\\end{cases}\n",
    "\\Leftrightarrow\n",
    "\\begin{cases}\n",
    "-2x_1 - x_2 = -2 \\\\\\\\\n",
    "4x_1 - x_2 = -8 \\\\\\\\\n",
    "-1x_1 - x_2 = -2\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Let's see their graphical representation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQkAAAECCAYAAADtryKnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHAxJREFUeJzt3XuUlPWd5/F3Q9MNmKYR0q2iIhjDt5HFGNozROMN1ImH\nHCaTM86ZzBidGLxg3Nn1kiySmeweHc24k6w7EhMnG5Q4juMGTpLNKJExRFFpvKwXTMgwX8ARdbm1\nZgYauWgjvX9UFamU1XV56rlV1ed1Dsfueqrr+VoU365P/Z7n+bYMDQ0hIjKcEUkXICLppiYhIiWp\nSYhISWoSIlKSmoSIlKQmISIltVZyJzN7Cdid/fY1d1+Qt+0q4GpgELjd3VeGXqWIJKal3HESZtYO\nrHP33iLbjgF+BswCxgJrgV53H4ygVhFJQCXvJD4GHGVm/wSMBP7c3Z/LbvsdYK27HwIGzGwzcBrw\nYiTVikjsKvlMYj/wDXf/FHAt8KCZ5X5uHLAn777vAJ3hligiSaqkSWwCHgRw983Ar4HjstsGyDSK\nnA5+89mFiDSASuLGF4GZwHVmNolMI9iR3fY8cJuZtQFjgB5gQ6kHGxoaGmppaQlecci2vLmbm5Y8\nxcTO0dz95TmMHT0q6ZJCs31gJ9c/egu9k2ay6JwvJV2OJC/QP7xKmsS9wDIzexo4TKZp/Gcz2+zu\nj5jZEjIfWLYAX3X390pW2dLCW2/tDVJrJDpHj+SSuR9l+epN3LNiPZdf3JN0SQB0dXXU/Dytfu1Z\nAGZ0nhracx5GXWFTTZXp6uoI9HNlm0R2peLzBTc/m7f9XjKNpG597qJp9L2yjTXrt9Pb082MKROS\nLikUL7/1C1pbRnJa16lJlyJ1TAdTAaNaR7Lg09MZ0dLC93+6kQPvHkq6pJrt2tfPtnd2MH3iNMa0\njkm6HKljahJZU44dx7wzJ/PrgXdZ8cSWpMup2Uv9vwTg412nJVyJ1Ds1iTzzz5rK8V1HsWb9dn61\n9d+SLqcmihoSFjWJPKNaRzRE7FDUkDCpSRRohNihqCFhUpMoot5jh6KGhElNooh6jh2KGhI2NYlh\n1GvsUNSQsKlJlFCPsUNRQ8KmJlFCvcUORQ2JgppEGfUUOxQ1JApqEhWol9ihqCFRUJOoQD3EDkUN\niYqaRIXSHjsUNSQqahJVSHPsUNSQqKhJVCGtsUNRQ6KkJlGlNMYORQ2JkppEAGmLHYoaEiU1iQDS\nFDsUNSRqahIBpSV2KGpI1CqdBdoNvABc6O6b8m6/AVgA9GdvuiY7m6MpzD9rKi9vfjvRC+gqakjU\nyr6TMLNW4G/JTPIqNAu4zN3nZv80TYOA5GOHoobEoZK48U3gHmB7kW29wGIze9rMbg61sjqRZOxQ\n1JA4lGwSZvYFoN/df0bx6T8PAQuBOcDZZjYv9ArrQFKrHYoaEody7ySuAC4ysyeA04G/y34+kXOX\nu/9bdqr4SuDjEdWZaknEDkUNiUvJDy7d/bzc19lGcY2792e/HwdsMLMe4AAwlwoneQUdNxalWmvq\n6urgkgsGWL56Ew8/+wbXXfKxSGt6qv9pAM77yOzYn89G/PuLQhprCqKi1Y2sIQAz+2PgKHdfamaL\ngTXAQeDn7r6qkgdK44zEMGq64PRJ9L2yjVXPbGXGSeNrWu0oV9ParS/Q2jKSKe1TY30+0zrjUjWV\nF9ks0Bx3n5v9clPebQ8CDwbacwPKxY7b7n+R7/90I7cumM2Y9mr6cGVyUWPmh6crakjkdDBVyOJY\n7dCqhsRJTSICUa92aFVD4qQmEYEoVzu0qiFxU5OISFSxQ1FD4qYmEaEoYoeihsRNTSJCYccORQ1J\ngppExMKMHYoakgQ1iRiEFTsUNSQJahIxCCN2KGpIUtQkYlJr7FDUkKSoScSoltihqCFJUZOIUdDY\noaghSVKTiFmQ2KGoIUlSk0hAtbFDUUOSpCaRgGpih6KGJE1NIiGVxg5FDUmamkSCKokdihqSNDWJ\nBJWLHYoakgZqEgkrFTsUNSQN1CRSYLjYoaghaVBRkzCzbjN7w8ymFdw+38yeN7M+M7symhIbX7HY\nsX1gp6KGpELgWaDZ2+8ELgTOB64uGNwjVSiMHc+8+RKgqCHJq2UW6HRgs7sPuPsgsBY4J+T6mkp+\n7Hhiy/OKGpIKtcwCHQfsyft+L9BZboe7fv54lSU2j1zsGDlmH/0HdzFt/EcVNSRx5SbHXAEcNrOL\n+M0s0N/LjvobINMocjqA3eV2uGXJt+meez4nX3MVI0ePDlh2+NIykq2rq4OZrw7i78Lh3ZNSU1c+\n1VSZNNYUROBZoMBG4BQzG0/m84pzgW+U2+GHTvkI/Y+vYffGTRy38Drajz8+ePUhSdtItr2tW+Hg\nCF5+sYU1U1+vaVxg2NL2XIFqqlTQplXNEuiRWaBmdmV2kviNwGNAH7DU3XeUe5CZd9zO+Asv4r0d\n23nj9lvY0/d0oMIb1a59/Wzft5OeCcaIw22xTSkXGU6ts0BXAiur2eGIUaPo/tyljJnWw65lS9m1\n7F4O+L/QfenljGhvr+ahGlLuAKoLp83m+L3tPLLudVY8sYXLL+5JuDJpVokdTNUxq5eT/uuttE+Z\nysC6Pt647Rbe3bYtqXJSI3cA1RmTTot8XKBIJRI94nJUVxcnLvqq4kdW/rkaY9vGRDouUKRSiR+W\nnYsfx33pz2gZOZJdy+5l533f4/C77yZdWuyKnasRx5RykVISbxI5ih/Dn6uh2CFJSk2TgOaOH6VO\nC1fskCSlqklA88aPcqeFK3ZIUlLXJHKaLX5Uclq4YockIbVNAponflR6BSrFDklCqpsENEf8qOYK\nVIodErfUN4mcRo4f1V6BSrFD4lQ3TQIaM34EuditYofEqa6aBDRe/Ah6sVvFDolL3TWJnEaJH7Vc\n7FaxQ+JQt00C6j9+1DpXQ7FD4lDXTQLqO36EMVdDsUOiVvdNIqce40dYczUUOyRKDdMkoL7iR5gj\n/BQ7JEoN1SSgfuJH2CP8FDskKg3XJHLSHj+iGOGn2CFRaNgmAemNH1FNC1fskCiUvRCumY0AvgcY\ncBhY6O7/nLf9BmABkLvU/jXuvjmCWgNJ44V3o5wWnosduoCuhKWSq2XPB4bc/WwzOw/4OvD7edtn\nAZe5+8tRFBiWjlm9jD5xMtu/+x0G1vVx8LXXEpv7EfW08PlnTeXlzW+zZv12enu6UzW3Q+pP2bjh\n7j8Brs5+OwX494K79AKLzexpM7s53PLClYb4EVXUyKfYIWGq6DMJdz9sZt8H7gIeLNj8ELAQmAOc\nbWbzQq0wZMVWPzbf9a3YVj+ijBr5tNohYWkZGhqq+M5m1g08D0x39wPZ28a5+0D262uBCe5+e4mH\nqXyHETu4axf+1/+Dd7a8ypgTTqBn0U2MnTw50n1+ZdVtbNu7i6Wf+WvGtkU7DHjw0Pvc8D+f5PWd\ne/nLa87k9Gndke5PUq9w6HdlP1SuSZjZ54ET3P0OMxsHvEymSbyX/X4D0AMcAJYD97r7qhIPOZSm\nGYmHBwfZ99P/w46HV9LS1kb3pZfR+clzItnXrn393PrcN5n54eksPO2KkvcNa5bk1p0D3Hb/ixzd\n0catC2Yzpr3ioW2R1hUm1VSZrq6OQE2ikrjxI+DjZvYk8ChwPfAH2XmgA8BiYA3wJLChTINInRGj\nRnHylV+M5eCruKJGPsUOqVXZXyvuvh/4oxLbH+SDn1PUnThWP6Je1RiOVjukFg19MFW1olz9iGNV\nYzha7ZBaqEkUiOrcjySiRj7FDglKTWIYYZ/7kVTUyKdzOyQINYkSwoofSUaNfIodEoSaRBlhxI+k\no0Y+xQ6plppEhWqJH2mIGvkUO6QaahJVCBI/0hI18il2SDXUJKpUbfxIU9TIp9ghlVKTCKjS+JG2\nqJFPsUMqoSZRg3LxI41RI59ih1RCTaJGpeJHWqNGPsUOKae2UwLliGLnfrz6yQ5a29MZNfLp3A4p\nRe8kQlQYP+b82Jmza1wqo0Y+xQ4pRU0iZLn4sfMPzuPwiBZOXe2pnPtRSLFDhqMmEZG1E3fzg3kf\nZtRJJ6Vy7kcxWu2QYtQkIpBb1TjhpFM56ea/SN3cj+EodkgxahIRyF/VqJexgzmKHVJITSICxQ6g\nSvvYwXyKHZJPTSJkpQ6gSsPcj0oodkg+NYmQlTuAql7ih2KH5IQxC3Q+8DVgEFjm7ksjqrUuVHqu\nRprGDg5HB1kJVPZO4sgsUDLN4Ou5DWbWCtwJXAicD1ydHeDTlKo9VyPt8UOxQ6D2WaDTgc3uPuDu\ng8BaIJrJNnUgyLkaw8WP9w8ejKrMqih2SK2zQMcBe/K+3wt0hlZdnanltPDC1Y9XblqUmtWP/NWO\n9Zv6ky5HYlbTLFAzmwnc4e6fzm6/E1jr7j8q8TCpmQUapu0DO7n+0VvonTSTRed8KfDjHB4cZOv9\nD7Dj4ZWMaGvj5IVXccwFc0OsNJgtb+7mpiVPMbFzNHd/eQ5jR49KuiSpXqAxf5V8cHlkFihwEHg/\n+wdgI3CKmY0H9gPnAt8o95gpnJFYc02rX3sWgBmdp9b8WB2f+UM6Z8xg013fYsuSb9P/wnq6L72c\nEe3tNT1uLTpHj2TeJybzyLrXuWfFei6/uCexWgqldO5mKmsKotZZoIeAG4HHgD5gqbvvCFRJnQv7\nClQTz5yduoOv5p81lZOO7dBBVk2mqrgRklRNFYfau34108Krrenw4CBv/3A5u1f/LPKp55XYc/B9\nbrrrqdCmlIchrb+1U1hTZFPFpYwor0CVtoOvTjlxvFY7moyaRAjiuNhtms790LkdzUVNokZxXuw2\nLQdf6SCr5qImUaO4L3ablvihg6yah5pEjZKaq5GG+KHY0RzUJGqQ9FyNpOOHYkdzUJOoQRrmaiQd\nPxQ7Gp+aRA3SNMIvyfih2NHY1CQCSjpqFJNU/FDsaGxqEgGlIWoUk1T8UOxoXGoSAaUpahSTRPxQ\n7GhMahIBpDFqFBN3/FDsaExqEgGkNWoUE3f8UOxoPGoSAaQ9ahQTZ/xQ7GgsahJVqpeoUUxc8UOx\no7GoSVSpnqJGMXHFD8WOxqEmUaV6jBrFxBE/FDsag5pEFeo5ahQTdfxQ7GgMahJVqPeoUUzU8UOx\no/6pSVShUaJGMVHGD8WO+laySZhZq5n9nZk9ZWbPZud+5m+/wcw2mNnj2T8fjbbc5DRa1Cgmqvih\n2FHfyr2T+DzwtrufC8wD7i7YPgu4zN3nZv9sjqLINGjEqFFMVPFDsaN+lWsSy8kMCYbM9J/Bgu29\nwGIze9rMbg67uDRp5KhRTBTxQ7GjPpVsEu6+3933mVkHsAL484K7PAQsBOYAZ5vZvGjKTFYzRI1i\nisWPXT9/PPjjKXbUpbLDeczsRDJTvO529/sLto1z94Hs19cCE9z99jL7rLtZoD/81U/5wYaH+Y+z\nv8C5U2YnXU4ifv3Mc2z+1t28v28/3XPP5+RrrmLk6NGBHuuBRzeyfPUmLj5zCtdd8rFwC5VSwp8F\nambHAP8EXOfuTxRsGwdsMLMe4AAwF7i3kp2mcLJRyZrWbn2B1paRTGmfGlvtqZsAdcqpTP6LW3jr\n3r+l//E17N64ieMWXkf78cdX/VAXnD6Jvle2seqZrcw4aTwzpkyoqbTUPVekt6Ygyn0msRgYD3zN\nzJ7IrmD8cXYO6EB2+xrgSWCDu68KVEWKNWvUKGZUVxcz77i95tUPxY76UvKdhLtfT2ZA8HDbHwQe\nDLuoNGmWVY1K5VY/xkzrYdeypexadi8H/F+qnnqeW+14ZN3rrHhiS6qmlMtv08FUZTTbqkalwlj9\n0GpHfVCTKEFRo7RaD75S7KgPahIlKGqUV+vBVzrIKv3UJEpQ1KhcLfFDsSPd1CSGoahRvaDxQ7Ej\n3dQkhqGoEUzQ+KHYkV5qEsNQ1KhNkPih2JFOahJFKGqEo9r4odiRTmoSRShqhKfa+KHYkT5qEkUo\naoSvmvih2JEuahIFFDWiU2n8UOxIFzWJAooa0ao0fih2pIeaRAFFjXhUEj8UO9JBTSKPoka8ysUP\nxY50UJPIo6gRv3LxQ7EjeWoSeRQ1klMqfih2JEtNIktRI3nDxQ/FjmSpSWQpaqTDcPFj8tHtih0J\nUZPIUtRIl2Lx41NT2hQ7EqAmAWwf2KmokUKF8WP7HX/JF44bUOyIWblL6rcC9wFTgDbgdnd/OG/7\nfDITvgaBZe6+NLpSo/PMmy8BihppVHjh3fd/9CBXfeQ07ts9QxfQjUngWaDZBnIncCFwPnC1mXVH\nVGeknn3zJUWNlMuPH0e/+gsW7FjFhuf+WbEjBrXMAp0ObHb3AXcfBNYC54RfYrR27evn9T3bFDXq\nQH78GH/g3/nT/7eSvgd+wv6DhSNqJUy1zAIdB+zJ+34v0Bl+idHSqkZ9yV/9GNHayvmvP8nPF/9V\nzVPPZXglP5OAD8wC/UHepgEyjSKnA9hdyU6DjhuLwi9f3EDriFbm9sxmbFu63kmk6XnKl4a6uj51\nPt3/wXj8v9zKxH99hVdv+W+c/rVFjJ08OenSjkjD8xSGwLNAgY3AKWY2HtgPnAt8o5KdpmVGYi5q\n9E6ayb49h9hHOuqCdM6ShJTV1fYhjv3yIp75m6Wc0b+R9Tcuovvzl9H5yeRTb6qep6wkZoEeAm4E\nHgP6gKXuviNQFQnJRY0zT+xNuBIJauoJE5nwuUv50bHnM0hL1XM/pLxaZ4GuBFaGXVRccgdQnTHp\nNPbt0Zp7vfrcRdPoe2Ub32s/mmvee4GBdX0cfO21wFPP5bc17cFU+edqpO2zCKnOqNaRLPj0dPa2\njWPZMRfxoTkX1jT1XH5b0zYJrWo0ltwp5W+9c4jVE88IPHZQPqhpm4TO1Wg8+aeUvzFhas1TzyWj\nKZuETgtvTIWnlB8ad3RNU88loymbhKJG4yq8klWtU8+lSZuEokZjK3Ylq1qmnje7pmsSihqNb7gr\nWQWdet7smq5JKGo0h+EuoKv4Ub2maxKKGs2j1AV0FT8q11RNQlGjuZS7gK7iR2WaqkkoajSfcnM7\nFD/Ka6omoajRnCqZ26H4MbymaRKKGs2r0rkdih/FNU2TUNRobpWOC1T8+KCmaRKKGlLNuEDFj99o\niiahqCFQ/ZRyxY+MpmgSihqSU+2UcsWPJmkSihqSL8iU8maOHw3fJBQ1pFDQKeXNGj8avkkoakgx\n1caOnGaMH2XnbgCY2WzgDnefU3D7DcACoD970zXuvjncEmujqCHDmX/WVF7e/DZr1m+nt6ebGVMm\nVPyzHbN6GX3iZLZ/9zsNf+Hdsu8kzOwrwPeA9iKbZwGXufvc7J9UNQhFDSklaOw48vNNEj8qiRtb\ngM8Os60XWGxmT5vZzeGVFQ5FDSknaOzIGS5+vH/wYATVJqNsk3D3HwPDtdiHgIXAHOBsM5sXYm01\nU9SQSgRZ7ShUuPrxyk2LGmb1o2VoaKjsnczsJOAhdz+r4PZx7j6Q/fpaYIK7317m4crvMATbB3Zy\n/aO30DtpJovO+VIcu5Q6tuXN3dy05Ckmdo7m7i/PYezoUYEe5/DgIFvvf4AdD69kRFsbJy+8imMu\nmBtytYG1BPmhij64LLYDMxsHbDCzHuAAMBe4t5IHimNG4urXngVgRuepZfeX1rmNaasJ0llXGDV1\njh7JvE9M5pF1r3PPivVcfnFP4Mfq+Mwf0jljBpvu+hZblnyb/hfW033p5YxoL/axXnyimgWabwgg\nbxboAJlZoWuAJ4EN7r4qUBURUNSQaoURO3Imnjm7YQ6+qihuhGwo6t9Eu/b1c+tz32Tmh6ez8LQr\nyt6/UX87RiGNdYVZ09adA9x2/4sc3dHGrQtmM6a9mjfbH6zp8OAgb/9wObtX/4yWtja6L01u6nlX\nV0eguNGQB1NpVUOCqnW1o1AjHHzVkE1CUUNqEWbsyKnncz8arknoACqpVa0HWQ37uHV68FXDNQlF\nDQlD2LEjpx7jR8M1CUUNCUsUsSOnnuJHQzUJRQ0JU1Sx48jj10n8aKgmoaghYYsqduTUQ/xoqCah\nqCFRiDJ25KQ5fjRMk1DUkKhEHTuO7Cel8aNhmoSihkQp6tiRk8b40TBNQlFDohZH7MhJU/xoiCah\nqCFxiCt2HNlfSuJHQzQJRQ2JS1yxIycN8aMhmoSihsQpztiRk2T8qPsmoaghcYs7dhzZb0Lxo+6b\nhKKGJCHu2JGTRPyo+yahqCFJSSJ25MQZP+q6SShqSJKSih1H9h9T/KjrJqGoIUlLKnbkxBE/KmoS\nZjbbzJ4ocvt8M3vezPrM7MrQqqqQooakQZKxIyfK+BF4zJ+ZtQJ3AhcC5wNXm1l3KFVVQFFD0iLp\n2HGkjojiRy1j/qYDm919wN0HgbVAbJcBVtSQNEk6duREET9qGfM3DtiT9/1eoDNwJVVS1JC0yY8d\n6zf1J1pLsfgRVC0fXA6QaRRH6gJ21/B4Fdu1/y22vbODngmKGpIe+bFjyfL1icWOI/UUxI+gqpkF\n+r/d/cy821qBXwGzgf3AOmC+u+8IXI2IpE4144mOjPkDjnL3pWZ2I/AYmTmhS9UgRBpPEmP+RKSO\n1PXBVCISPTUJESlJTUJESlKTEJGSqlndCMTMxgF/T+aYilHATe7+bMF9rgKuBgaB2919ZdR1Zff7\nWeASd7+0yLa7gLPIHCQG8Bl331t4v5hrivV5MrPRZP7uuskcF/On7v7rgvv8BJiQremAu386olpa\ngO8AHwMOAle6+7/mbU/qNVSurqReR7OBO9x9TsHt84GvkXmelrn70nKPFXmTAG4EVrv7EjObBjwE\n9OY2mtkxwJ8Bs4CxwFozeyx7qHdkzOxvgN8F1g9zl1nAp9w9tjN2StWU0PN0LfALd7/VzP6IzIvr\n+oL7nOLuMyKsIef3gXZ3Pyv7D+DO7G2JvYbK1ZWVxOvoK8BlwDsFt+fOt+oFDgB9ZvaP7l7y8NA4\n4sadwHezX48iU1y+3wHWuvshdx8ANgNxnJDRR+YfwQdkfzt8FPhfZrbWzK6IoZ6SNZHM83Q2sCr7\n9aNkTuY7IntC33gz+0cze8rMInkXUViLuz8HnJG3LanXUMm6EnwdhXq+VajvJMzsi8ANZA68asn+\n9wp3f9HMjgUeAP5TwY8VngPyDiGeA1KiphVmdt4wP3YUsIRMg2sFnjCz/+vuGxKsKc7niWxdO/P2\nuZffPgwfoA34JnAXMJHMb6bn3P3tsOrKU/j/f8jMRrj74SLbQn1uaqgr0tfRcNz9x9mjpMvVWtH5\nVqE2CXe/D7iv8HYzmwn8A5nPI9YWbI70HJDhaipjP7DE3Q8CmNnjZDJnKH+5AWuK/Xkysx9m9zPc\n/nYC383+g3jLzF4GDIiiSQzk1QKQ+4eY25bIeURl6or0dRRAoOcp8rhhZqcCy4E/cffHitzleeBs\nM2szs06gh+SexJxpZHJti5mNIvOW8qWEa0rieeoD5mW/ngcUXpzgQjJ/t5jZh4AZwMaoazGzTwC/\nzNuW5GuoVF1Jv45aCr7fCJxiZuPNrA04F3im3IPE8cHl18lcsOaubEbb7e6fNbMbyOSjR8xsCZl8\n1AJ81d3fi6GuDyio6e+B54D3gPvdPaoXfzU1xf083QPcb2ZPA+8Cf5Kt6b8DK9x9lZn9rpk9A7wP\nLI7wA7ofAxeZWV/2+ytS8hoqV1eSr6NQzrfSuRsiUpIOphKRktQkRKQkNQkRKUlNQkRKUpMQkZLU\nJESkJDUJESlJTUJESvr/JWqSL6j9kw4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108e39b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x1 = np.linspace(-5, 5, 1000)\n",
    "x2_1 = -2*x1 + 2\n",
    "x2_2 = 4*x1 + 8\n",
    "x2_3 = -1*x1 + 2\n",
    "\n",
    "plt.plot(x1, x2_1)\n",
    "plt.plot(x1, x2_2)\n",
    "plt.plot(x1, x2_3)\n",
    "plt.xlim(-2., 1)\n",
    "plt.ylim(1, 5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We actually see that there is no solution.\n",
    "\n",
    "Putting this into the matrix form we have:\n",
    "\n",
    "$$\n",
    "\\bs{A}=\n",
    "\\begin{bmatrix}\n",
    "    -2 & -1 \\\\\\\\\n",
    "    4 & -1 \\\\\\\\\n",
    "    -1 & -1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\bs{x}=\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\\\\\n",
    "    x_2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\bs{b}=\n",
    "\\begin{bmatrix}\n",
    "    -2 \\\\\\\\\n",
    "    -8 \\\\\\\\\n",
    "    -2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "So we have:\n",
    "\n",
    "$$\n",
    "\\bs{Ax} = \\bs{b}\n",
    "\\Leftrightarrow\n",
    "\\begin{bmatrix}\n",
    "    -2 & -1 \\\\\\\\\n",
    "    4 & -1 \\\\\\\\\n",
    "    -1 & -1\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\\\\\n",
    "    x_2\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "    -2 \\\\\\\\\n",
    "    -8 \\\\\\\\\n",
    "    -2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "We will now calculate the pseudoinverse of $\\bs{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.11290323,  0.17741935, -0.06451613],\n",
       "       [-0.37096774, -0.27419355, -0.35483871]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[-2, -1], [4, -1], [-1, -1]])\n",
    "A_plus = np.linalg.pinv(A)\n",
    "A_plus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have calculated the pseudoinverse of $\\bs{A}$:\n",
    "\n",
    "$$\n",
    "\\bs{A}^+=\n",
    "\\begin{bmatrix}\n",
    "    -0.11290323 &  0.17741935 & -0.06451613 \\\\\\\\\n",
    "    -0.37096774 & -0.27419355 & -0.35483871\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "we can use it to find $\\bs{x}$ knowing that:\n",
    "\n",
    "$$\n",
    "\\bs{x}=\\bs{A}^+\\bs{b}\n",
    "$$\n",
    "\n",
    "with:\n",
    "\n",
    "$$\n",
    "\\bs{x}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "    x1 \\\\\\\\\n",
    "    x2\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.06451613],\n",
       "       [ 3.64516129]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.array([[-2], [-8], [-2]])\n",
    "res = A_plus.dot(b)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\bs{A}^+\\bs{b}&=\n",
    "\\begin{bmatrix}\n",
    "    -0.11290323 &  0.17741935 & -0.06451613 \\\\\\\\\n",
    "    -0.37096774 & -0.27419355 & -0.35483871\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "    -2 \\\\\\\\\n",
    "    -8 \\\\\\\\\n",
    "    -2\n",
    "\\end{bmatrix}\\\\\\\\\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "    -1.06451613 \\\\\\\\\n",
    "    3.64516129\n",
    "\\end{bmatrix}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "In our two dimensions, the coordinates of $\\bs{x}$ are\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    -1.06451613 \\\\\\\\\n",
    "    3.64516129\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Let's plot this point along with the equations lines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQkAAAECCAYAAADtryKnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHIVJREFUeJzt3XuUlPWd5/F3Q9MNGBoEu1VULsbwbWQxhvYM0agR1MRD\nDut4xjmTGaMTgxeMO7tekkUzye7R0aybZN2RmDgZUOI4rhs9MZtRImNQUcHbqmBChvkCRtDl1poJ\nNHLRRnr/qCpSltV1eeq5VdXndQ7H7nqqn/paFN+uT/2e5/m2DAwMICIymCFJFyAi6aYmISIlqUmI\nSElqEiJSkpqEiJSkJiEiJbVWciczexXYmf32DXefl7ftcuAKoB+41d2Xhl6liCSmpdxxEmbWDjzn\n7j1Fth0J/BKYAYwEVgI97t4fQa0ikoBK3kl8EjjMzP4ZGAr8tbu/mN32R8BKdz8A9JnZBuAk4JVI\nqhWR2FXymcRe4Lvu/nngKuB+M8v9XAewK+++7wKjwy1RRJJUSZNYD9wP4O4bgN8BR2e39ZFpFDmj\n+MNnFyLSACqJG18BpgNXm9l4Mo1gW3bbS8AtZtYGjAC6gbWldjYwMDDQ0tISvOKQbXxrJ9cvfIZx\no4dz59dmMXL4sKRLCs3Wvu1c89hN9IyfzoIzvpp0OZK8QP/wKmkSdwNLzOxZ4CCZpvGfzGyDuz9q\nZgvJfGDZAnzD3d8vWWVLC2+/vTtIrZEYPXwoF87+BA8uX89dD63hkvO6ky4JgM7OUTU/T8vfeAGA\naaNPDO05D6OusKmmynR2jgr0c2WbRHal4ksFN7+Qt/1uMo2kbn3x3Cmsem0LK9Zspae7i2mTxiZd\nUihWv/0rWluGclLniUmXInVMB1MBw1qHMu8LUxnS0sKPf7GOfe8dSLqkmu3Y08uWd7cxddwURrSO\nSLocqWNqElmTjupgzqkT+F3fezz01Maky6nZq72/BuBTnSclXInUOzWJPHNPm8wxnYexYs1WfrPp\n35IupyaKGhIWNYk8w1qHNETsUNSQMKlJFGiE2KGoIWFSkyii3mOHooaESU2iiHqOHYoaEjY1iUHU\na+xQ1JCwqUmUUI+xQ1FDwqYmUUK9xQ5FDYmCmkQZ9RQ7FDUkCmoSFaiX2KGoIVFQk6hAPcQORQ2J\nippEhdIeOxQ1JCpqElVIc+xQ1JCoqElUIa2xQ1FDoqQmUaU0xg5FDYmSmkQAaYsdihoSJTWJANIU\nOxQ1JGpqEgGlJXYoakjUKp0F2gW8DJzj7uvzbr8WmAf0Zm+6MjuboynMPW0yqze8k+gFdBU1JGpl\n30mYWSvwd2QmeRWaAVzs7rOzf5qmQUDysUNRQ+JQSdz4HnAXsLXIth7gRjN71sxuCLWyOpFk7FDU\nkDiUbBJm9mWg191/SfHpPw8A84FZwOlmNif0CutAUqsdihoSh3LvJC4FzjWzp4CTgX/Ifj6Rc4e7\n/1t2qvhS4FMR1ZlqScQORQ2JS8kPLt39s7mvs43iSnfvzX7fAaw1s25gHzCbCid5BR03FqVaa+rs\nHMWFZ/fx4PL1PPLCm1x94ScjremZ3mcB+OzHZ8b+fDbi318U0lhTEBWtbmQNAJjZnwOHuftiM7sR\nWAHsB55w92WV7CiNMxLDqOnsk8ez6rUtLHt+E9MmjqlptaNcTSs3vUxry1AmtU+O9flM64xL1VRe\nZLNAc9x9dvbL9Xm33Q/cH+iRG1Audtxy7yv8+BfruHneTEa0V9OHK5OLGtOPmKqoIZHTwVQhi2O1\nQ6saEic1iQhEvdqhVQ2Jk5pEBKJc7dCqhsRNTSIiUcUORQ2Jm5pEhKKIHYoaEjc1iQiFHTsUNSQJ\nahIRCzN2KGpIEtQkYhBW7FDUkCSoScQgjNihqCFJUZOISa2xQ1FDkqImEaNaYoeihiRFTSJGQWOH\nooYkSU0iZkFih6KGJElNIgHVxg5FDUmSmkQCqokdihqSNDWJhFQaOxQ1JGlqEgmqJHYoakjS1CQS\nVC52KGpIGqhJJKxU7FDUkDRQk0iBwWKHooakQUVNwsy6zOxNM5tScPtcM3vJzFaZ2WXRlNj4isWO\nrX3bFTUkFQLPAs3efjtwDnAWcEXB4B6pQmHseP6tVwFFDUleLbNApwIb3L3P3fuBlcAZIdfXVPJj\nx1MbX1LUkFSoZRZoB7Ar7/vdwOhyD7jjiSerLLF55GLH0BF76N2/gyljPqGoIYkrNznmUuCgmZ3L\nH2aB/vvsqL8+Mo0iZxSws9wDblz4A7pmn8XxV17O0OHDA5YdvrSMZOvsHMX01/vx9+DgzvGpqSuf\naqpMGmsKIvAsUGAdcIKZjSHzecWZwHfLPeDHTvg4vU+uYOe69Rw9/2rajzkmePUhSdtItt2tm2D/\nEFa/0sKKyZtrGhcYtrQ9V6CaKhW0aVWzBHpoFqiZXZadJH4d8DiwCljs7tvK7WT6bbcy5pxzeX/b\nVt689SZ2rXo2UOGNaseeXrbu2U73WGPIwbbYppSLDKbWWaBLgaXVPOCQYcPo+uJFjJjSzY4li9mx\n5G72+b/SddElDGlvr2ZXDSl3ANU5U2ZyzO52Hn1uMw89tZFLzutOuDJpVokdTDVqRg8T/8vNtE+a\nTN9zq3jzlpt4b8uWpMpJjdwBVKeMPynycYEilUj0iMthnZ0ct+Abih9Z+edqjGwbEem4QJFKJX5Y\ndi5+HP3Vv6Jl6FB2LLmb7fcs4uB77yVdWuyKnasRx5RykVISbxI5ih+Dn6uh2CFJSk2TgOaOH6VO\nC1fskCSlqklA88aPcqeFK3ZIUlLXJHKaLX5Uclq4YockIbVNAponflR6BSrFDklCqpsENEf8qOYK\nVIodErfUN4mcRo4f1V6BSrFD4lQ3TQIaM34EuditYofEqa6aBDRe/Ah6sVvFDolL3TWJnEaJH7Vc\n7FaxQ+JQt00C6j9+1DpXQ7FD4lDXTQLqO37kosbqR/r53Oee4PLLH+b3vy97ca8PUeyQqFV8PYm0\nGzWjh+HHTWDrj35I33Or2P/GG6m58tVgVr/9KwY+aOGRv7+KA/tGsmbNAHAfixZdUNV+5p42mdUb\n3mHFmq30dHel6kpWUv/q/p1EvnqKH7mosee3R3Bg38jsrS1s3txR8ueKUeyQKDVUk4D6iR+5qDH8\nnQ6yVwYEBpg4sS/Q/hQ7JCoNEzcKpT1+5FY1vnnl+RzcdB+bN3cwcWIf3/nOrMD7VOyQKDTcO4l8\naY0f+asaRx9xFIsWXcDjj5/NokUXcPjhYwLvV7FDolD2nYSZDQEWAQYcBOa7+7/kbb8WmAfkLrV/\npbtviKDWQNJ44d0op4XnYocuoCthqSRuzAUG3P10M/ss8G3gj/O2zwAudvfVURQYljTFj6inhSt2\nSJjKxg13/zlwRfbbScDvC+7SA9xoZs+a2Q3hlheuNMSPWg+gqoRih4Spos8k3P2gmf0YuAO4v2Dz\nA8B8YBZwupnNCbXCkBVb/dhwx/djW/2IMmrk02qHhKVlYGCg/L2yzKwLeAmY6u77srd1uHtf9uur\ngLHufmuJ3VT+gBHbv2MH/p3/wbsbX2fEscfSveB6Rk6YEOljfn3ZLWzZvYPF53+HkW3RDgPuP/AB\n1/7Pp9m8fTd/c+WpnDylK9LHk9QrHPpd2Q+VaxJm9iXgWHe/zcw6gNVkmsT72e/XAt3APuBB4G53\nX1ZilwNpmpF4sL+fPb/4P2x7ZCktbW10XXQxoz9zRiSPtWNPLze/+D2mHzGV+SddWvK+Yc2S3LS9\nj1vufYXDR7Vx87yZjGivbdU7rTMuVVN5nZ2jAjWJSuLGw8CnzOxp4DHgGuBPsvNA+4AbgRXA08Da\nMg0idYYMG8bxl30lloOv4ooa+RQ7pFZlf624+17gz0psv5+Pfk5Rd+JY/Yh6VWMwWu2QWjT0wVTV\ninL1I45VjcFotUNqoSZRIKpzP5KIGvkUOyQoNYlBhH3lq6SiRj5dyUqCUJMoIaz4kWTUyKfYIUGo\nSZQRRvxIOmrkU+yQaqlJVKiW+JGGqJFPsUOqoSZRhSDxIy1RI59ih1RDTaJK1caPNEWNfIodUik1\niYAqjR9pixr5FDukEmoSNSgXP9IYNfIpdkgl1CRqVCp+pDVq5FPskHIa9kK4cSt27sfrnxlFa3s6\no0Y+ndshpeidRIgK48esnzmzdnSkMmrkU+yQUtQkQpaLH9v/5LMcHNLCics9lXM/Cil2yGDUJCKy\nctxOfjLnCIZNnFg3U8+12iHFqElEILeqcezEE5l4wzdTN/djMIodUoyaRATyVzXqZexgjmKHFFKT\niECxA6jCPvU8Soodkk9NImSlDqBKw9yPSih2SD41iZCVO4CqXuKHYofkhDELdC7wLaAfWOLuiyOq\ntS5Ueq5GmsYODkYHWQlU9k7i0CxQMs3g27kNZtYK3A6cA5wFXJEd4NOUqj1XI+3xQ7FDoPZZoFOB\nDe7e5+79wEogmsk2dSDIuRqDxY8P9u+PqsyqKHZIrbNAO4Bded/vBkaHVl2dqeW08MLVj9euX5Ca\n1Y/81Y4163uTLkdiVtMsUDObDtzm7l/Ibr8dWOnuD5fYTWpmgYZpa992rnnsJnrGT2fBGV8NvJ+D\n/f1suvc+tj2ylCFtbRw//3KOPHt2iJUGs/GtnVy/8BnGjR7OnV+bxcjhw5IuSaoXaMxfJR9cHpoF\nCuwHPsj+AVgHnGBmY4C9wJnAd8vtM4UzEmuuafkbLwAwbfSJNe9r1Pl/yuhp01h/x/fZuPAH9L68\nhq6LLmFIe3tN+63F6OFDmfPpCTz63GbuemgNl5zXnVgthVI6dzOVNQVR6yzQA8B1wOPAKmCxu28L\nVEmdC/sKVONOnZm6g6/mnjaZiUeN0kFWTaaquBGSVE0Vh9q7fjXTwqut6WB/P+/89EF2Lv9l5FPP\nK7Fr/wdcf8czoU0pD0Naf2unsKbIpopLGVFegSptB1+dcNwYrXY0GTWJEMRxsds0nfuhczuai5pE\njeK82G1aDr7SQVbNRU2iRnFf7DYt8UMHWTUPNYkaJTVXIw3xQ7GjOahJ1CDpuRpJxw/FjuagJlGD\nNMzVSDp+KHY0PjWJGqRphF+S8UOxo7GpSQSUdNQoJqn4odjR2NQkAkpD1Cgmqfih2NG41CQCSlPU\nKCaJ+KHY0ZjUJAJIY9QoJu74odjRmNQkAkhr1Cgm7vih2NF41CQCSHvUKCbO+KHY0VjUJKpUL1Gj\nmLjih2JHY1GTqFI9RY1i4oofih2NQ02iSvUYNYqJI34odjQGNYkq1HPUKCbq+KHY0RjUJKpQ71Gj\nmKjjh2JH/VOTqEKjRI1ioowfih31rWSTMLNWM/sHM3vGzF7Izv3M336tma01syezfz4RbbnJabSo\nUUxU8UOxo76VeyfxJeAddz8TmAPcWbB9BnCxu8/O/tkQRZFp0IhRo5io4odiR/0q1yQeJDMkGDLT\nf/oLtvcAN5rZs2Z2Q9jFpUkjR41ioogfih31qWSTcPe97r7HzEYBDwF/XXCXB4D5wCzgdDObE02Z\nyWqGqFFMsfix44kng+9PsaMulR3OY2bHkZnidae731uwrcPd+7JfXwWMdfdbyzxm3c0C/elvfsFP\n1j7Cf5j5Zc6cNDPpchLxu+dfZMP37+SDPXvpmn0Wx195OUOHDw+0r/seW8eDy9dz3qmTuPrCT4Zb\nqJQS/ixQMzsS+Gfgand/qmBbB7DWzLqBfcBs4O5KHjSFk41K1rRy08u0tgxlUvvk2GpP3QSoE05k\nwjdv4u27/47eJ1ewc916jp5/Ne3HHFP1rs4+eTyrXtvCsuc3MW3iGKZNGltTaal7rkhvTUGU+0zi\nRmAM8C0zeyq7gvHn2TmgfdntK4CngbXuvixQFSnWrFGjmGGdnUy/7daaVz8UO+pLyXcS7n4NmQHB\ng22/H7g/7KLSpFlWNSqVW/0YMaWbHUsWs2PJ3ezzf6166nlutePR5zbz0FMbUzWlXD5MB1OV0Wyr\nGpUKY/VDqx31QU2iBEWN0mo9+Eqxoz6oSZSgqFFerQdf6SCr9FOTKEFRo3K1xA/FjnRTkxiEokb1\ngsYPxY50U5MYhKJGMEHjh2JHeqlJDEJRozZB4odiRzqpSRShqBGOauOHYkc6qUkUoagRnmrjh2JH\n+qhJFKGoEb5q4odiR7qoSRRQ1IhOpfFDsSNd1CQKKGpEq9L4odiRHmoSBRQ14lFJ/FDsSAc1iTyK\nGvEqFz8UO9JBTSKPokb8ysUPxY7kqUnkUdRITqn4odiRLDWJLEWN5A0WPxQ7kqUmkaWokQ6DxY8J\nh7crdiRETSJLUSNdisWPz09qU+xIgJoEsLVvu6JGChXGj623/Q1fPrpPsSNm5S6p3wrcA0wC2oBb\n3f2RvO1zyUz46geWuPvi6EqNzvNvvQooaqRR4YV3P3j4fi7/+Encs3OaLqAbk8CzQLMN5HbgHOAs\n4Aoz64qozki98Narihoplx8/Dn/9V8zbtoy1L/6LYkcMapkFOhXY4O597t4PrATOCL/EaO3Y08vm\nXVsUNepAfvwYs+/3/OX/W8qq+37O3v2FI2olTLXMAu0AduV9vxsYHX6J0dKqRn3JX/0Y0trKWZuf\n5okb/1vNU89lcCU/k4CPzAL9Sd6mPjKNImcUsLOSBw06biwKv35lLa1DWpndPZORbel6J5Gm5ylf\nGurq/PxZdP0748n/fDPjfvsar9/0Xzn5WwsYOWFC0qUdkobnKQyBZ4EC64ATzGwMsBc4E/huJQ+a\nlhmJuajRM346e3YdYA/pqAvSOUsSUlZX28c46msLeP5vF3NK7zrWXLeAri9dzOjPJJ96U/U8ZSUx\nC/QAcB3wOLAKWOzu2wJVkZBc1Dj1uJ6EK5GgJh87jrFfvIiHjzqLflqqnvsh5dU6C3QpsDTsouKS\nO4DqlPEnsWeX1tzr1RfPncKq17awqP1wrnz/ZfqeW8X+N94IPPVcPqxpD6bKP1cjbZ9FSHWGtQ5l\n3hemsrutgyVHnsvHZp1T09Rz+bCmbRJa1WgsuVPK3373AMvHnRJ47KB8VNM2CZ2r0XjyTyl/c+zk\nmqeeS0ZTNgmdFt6YCk8pP9BxeE1TzyWjKZuEokbjKrySVa1Tz6VJm4SiRmMrdiWrWqaeN7umaxKK\nGo1vsCtZBZ163uyarkkoajSHwS6gq/hRvaZrEooazaPUBXQVPyrXVE1CUaO5lLuAruJHZZqqSShq\nNJ9yczsUP8prqiahqNGcKpnbofgxuKZpEooazavSuR2KH8U1TZNQ1GhulY4LVPz4qKZpEooaUs24\nQMWPP2iKJqGoIVD9lHLFj4ymaBKKGpJT7ZRyxY8maRKKGpIvyJTyZo4fDd8kFDWkUNAp5c0aPxq+\nSShqSDHVxo6cZowfZeduAJjZTOA2d59VcPu1wDygN3vTle6+IdwSa6OoIYOZe9pkVm94hxVrttLT\n3cW0SWMr/tlRM3oYftwEtv7ohw1/4d2y7yTM7OvAIqC9yOYZwMXuPjv7J1UNQlFDSgkaOw79fJPE\nj0rixkbggkG29QA3mtmzZnZDeGWFQ1FDygkaO3IGix8f7N8fQbXJKNsk3P1nwGAt9gFgPjALON3M\n5oRYW80UNaQSQVY7ChWufrx2/YKGWf1oGRgYKHsnM5sIPODupxXc3uHufdmvrwLGuvutZXZX/gFD\nsLVvO9c8dhM946ez4IyvxvGQUsc2vrWT6xc+w7jRw7nza7MYOXxYoP0c7O9n0733se2RpQxpa+P4\n+Zdz5NmzQ642sJYgP1TRB5fFHsDMOoC1ZtYN7ANmA3dXsqM4ZiQuf+MFAKaNPrHs46V1bmPaaoJ0\n1hVGTaOHD2XOpyfw6HObueuhNVxyXnfgfY06/08ZPW0a6+/4PhsX/oDel9fQddElDGkv9rFefKKa\nBZpvACBvFmgfmVmhK4CngbXuvixQFRFQ1JBqhRE7csadOrNhDr6qKG6EbCDq30Q79vRy84vfY/oR\nU5l/0qVl79+ovx2jkMa6wqxp0/Y+brn3FQ4f1cbN82Yyor2aN9sfrelgfz/v/PRBdi7/JS1tbXRd\nlNzU887OUYHiRkMeTKVVDQmq1tWOQo1w8FVDNglFDalFmLEjp57P/Wi4JqEDqKRWtR5kNeh+6/Tg\nq4ZrEooaEoawY0dOPcaPhmsSihoSlihiR049xY+GahKKGhKmqGLHof3XSfxoqCahqCFhiyp25NRD\n/GioJqGoIVGIMnbkpDl+NEyTUNSQqEQdOw49TkrjR8M0CUUNiVLUsSMnjfGjYZqEooZELY7YkZOm\n+NEQTUJRQ+IQV+w49HgpiR8N0SQUNSQuccWOnDTEj4ZoEooaEqc4Y0dOkvGj7puEoobELe7Ycehx\nE4ofdd8kFDUkCXHHjpwk4kfdNwlFDUlKErEjJ874UddNQlFDkpRU7Dj0+DHFj7puEooakrSkYkdO\nHPGjoiZhZjPN7Kkit881s5fMbJWZXRZaVRVS1JA0SDJ25EQZPwKP+TOzVuB24BzgLOAKM+sKpaoK\nKGpIWiQdOw7VEVH8qGXM31Rgg7v3uXs/sBKI7TLAihqSJknHjpwo4kctY/46gF153+8GRgeupEqK\nGpI2+bFjzfreRGspFj+CquWDyz4yjeJQXcDOGvZXsR1732bLu9voHquoIemRHzsWPrgmsdhxqJ6C\n+BFUNbNA/7e7n5p3WyvwG2AmsBd4Dpjr7tsCVyMiqVPNeKJDY/6Aw9x9sZldBzxOZk7oYjUIkcaT\nxJg/EakjdX0wlYhET01CREpSkxCRktQkRKSkalY3AjGzDuAfyRxTMQy43t1fKLjP5cAVQD9wq7sv\njbqu7ONeAFzo7hcV2XYHcBqZg8QAznf33YX3i7mmWJ8nMxtO5u+ui8xxMX/p7r8ruM/PgbHZmva5\n+xciqqUF+CHwSWA/cJm7/zZve1KvoXJ1JfU6mgnc5u6zCm6fC3yLzPO0xN0Xl9tX5E0CuA5Y7u4L\nzWwK8ADQk9toZkcCfwXMAEYCK83s8eyh3pExs78FPgesGeQuM4DPu3tsZ+yUqimh5+kq4FfufrOZ\n/RmZF9c1Bfc5wd2nRVhDzh8D7e5+WvYfwO3Z2xJ7DZWrKyuJ19HXgYuBdwtuz51v1QPsA1aZ2T+5\ne8nDQ+OIG7cDP8p+PYxMcfn+CFjp7gfcvQ/YAMRxQsYqMv8IPiL72+ETwN+b2UozuzSGekrWRDLP\n0+nAsuzXj5E5me+Q7Al9Y8zsn8zsGTOL5F1EYS3u/iJwSt62pF5DJetK8HUU6vlWob6TMLOvANeS\nOfCqJfvfS939FTM7CrgP+I8FP1Z4Dsi7hHgOSImaHjKzzw7yY4cBC8k0uFbgKTP7v+6+NsGa4nye\nyNa1Pe8xd/Phw/AB2oDvAXcA48j8ZnrR3d8Jq648hf//B8xsiLsfLLIt1OemhroifR0Nxt1/lj1K\nulytFZ1vFWqTcPd7gHsKbzez6cD/IvN5xMqCzZGeAzJYTWXsBRa6+34AM3uSTOYM5S83YE2xP09m\n9tPs4wz2eNuBH2X/QbxtZqsBA6JoEn15tQDk/iHmtiVyHlGZuiJ9HQUQ6HmKPG6Y2YnAg8BfuPvj\nRe7yEnC6mbWZ2Wigm+SexJwpZHJti5kNI/OW8tWEa0rieVoFzMl+PQcovDjBOWT+bjGzjwHTgHVR\n12JmnwZ+nbctyddQqbqSfh21FHy/DjjBzMaYWRtwJvB8uZ3E8cHlt8lcsOaObEbb6e4XmNm1ZPLR\no2a2kEw+agG+4e7vx1DXRxTU9I/Ai8D7wL3uHtWLv5qa4n6e7gLuNbNngfeAv8jW9N+Bh9x9mZl9\nzsyeBz4AbozwA7qfAeea2ars95em5DVUrq4kX0ehnG+lczdEpCQdTCUiJalJiEhJahIiUpKahIiU\npCYhIiWpSYhISWoSIlKSmoSIlPT/AZJPuuHmDpkrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108fab2d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x1, x2_1)\n",
    "plt.plot(x1, x2_2)\n",
    "plt.plot(x1, x2_3)\n",
    "plt.xlim(-2., 1)\n",
    "plt.ylim(1, 5)\n",
    "\n",
    "plt.scatter(res[0], res[1])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe you would have expected the point being at the barycenter of the triangle (cf. [Least square solution in the triangle center](https://math.stackexchange.com/questions/471812/is-the-least-squares-solution-to-an-overdetermined-system-a-triangle-center)). This is not the case becase the equations are not scaled the same way. Actually the point is at the intersection of the three [symmedians](https://en.wikipedia.org/wiki/Symmedian) of the triangle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3.\n",
    "\n",
    "This method can also be used to fit a line to a set of points. Let's take the following data points:\n",
    "\n",
    "<img src=\"images/dataset-representation.png\" width=\"300\" alt=\"Representation of a set of data points\" title=\"Some data points\">\n",
    "<em>We want to fit a line to this set of data points</em>\n",
    "\n",
    "We have this set of $\\bs{x}$ and $\\bs{y}$ and we are looking for the line $y=mx+b$ that minimizes the error. The error can be evaluated as the sum of the differences between the fit and the actual data points. We can represent the data points with a matrix equations:\n",
    "\n",
    "$$\n",
    "\\bs{Ax} = \\bs{b}\n",
    "\\Leftrightarrow\n",
    "\\begin{bmatrix}\n",
    "    0 & 1 \\\\\\\\\n",
    "    1 & 1 \\\\\\\\\n",
    "    2 & 1 \\\\\\\\\n",
    "    3 & 1 \\\\\\\\\n",
    "    3 & 1 \\\\\\\\\n",
    "    4 & 1\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "    m \\\\\\\\\n",
    "    b\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "    2 \\\\\\\\\n",
    "    4 \\\\\\\\\n",
    "    0 \\\\\\\\\n",
    "    2 \\\\\\\\\n",
    "    5 \\\\\\\\\n",
    "    3\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Note that here the matrix $\\bs{A}$ represents the values of the coefficients. The column of 1 correspond to the intercepts (without it the fit would have the constraint to cross the origin). It gives the following set of equations:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "    0m + 1b = 2 \\\\\\\\\n",
    "    1m + 1b = 4 \\\\\\\\\n",
    "    2m + 1b = 0 \\\\\\\\\n",
    "    3m + 1b = 2 \\\\\\\\\n",
    "    3m + 1b = 5 \\\\\\\\\n",
    "    4m + 1b = 3\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "We have the set of equations $mx+b=y$. The ones are used to give back the intercept parameter. For instance, in the first equation corresponding to the first point we have well $x=0$ and $y=2$. This can be confusing because here the vector $\\bs{x}$ corresponds to the coefficients. This is because the problem is different from the other examples: we are looking for the coefficients of a line and not for $x$ and $y$ unknowns. We kept this notation to indicate the similarity with the last examples.\n",
    "\n",
    "So we will construct these matrices and try to use the pseudoinverse to find the equation of the line minimizing the error (difference between the line and the actual data points).\n",
    "\n",
    "Let's start with the creation of the matrix of $\\bs{A}$ and $\\bs{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1],\n",
       "       [1, 1],\n",
       "       [2, 1],\n",
       "       [3, 1],\n",
       "       [3, 1],\n",
       "       [4, 1]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[0, 1], [1, 1], [2, 1], [3, 1], [3, 1], [4, 1]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2],\n",
       "       [4],\n",
       "       [0],\n",
       "       [2],\n",
       "       [5],\n",
       "       [3]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.array([[2], [4], [0], [2], [5], [3]])\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the pseudoinverse of $\\bs{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -2.00000000e-01,  -1.07692308e-01,  -1.53846154e-02,\n",
       "          7.69230769e-02,   7.69230769e-02,   1.69230769e-01],\n",
       "       [  6.00000000e-01,   4.00000000e-01,   2.00000000e-01,\n",
       "          4.16333634e-17,   4.16333634e-17,  -2.00000000e-01]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_plus = np.linalg.pinv(A)\n",
    "A_plus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and apply it to the result to find the coefficients with the formula:\n",
    "\n",
    "$$\n",
    "\\bs{x}=\\bs{A}^+\\bs{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.21538462],\n",
       "       [ 2.2       ]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coefs = A_plus.dot(b)\n",
    "coefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the parameters of the fit. The slope is $m=0.21538462$ and the intercept is $b=2.2$. We will plot the data points and the regression line: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD9CAYAAAB3NXH8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEqdJREFUeJzt3X1s3Hdhx/G3nx/PvktyjvPsh9hftw2kUFKewroiEGMa\nDSD2IK1DUKYNRjWkSd2GEPwxCUQWNm0SHUgICpNSCZhKA1s1MqSIpQWBUEC0dP7avnPOSRPXjrnz\nQ/zs++2Pu/zipLbvF8f35O/nJVX19R6+37v67d/3fvfwq/A8DxFxR2WxJyAihaXoRRyj6EUco+hF\nHKPoRRyj6EUcU73VN7i8vOIlk7NbfbOBRSKNuDq+y/dd40M0GqoIcrkt39JXV1dt9U1q/DIYW+MX\nf/ygtLwXcYyiF3GMohdxjKIXcYyiF3GMohdxjKIXcYyiF3GMohdxjKIXcYyiF3GMohdxTKBP2Rlj\nLgCp7Mlha+3H8jclEcmnnNEbY+oAz1r7zgLMR8pUfyLJ6OQC7a11xZ6K5BBkeX8UaDLG/NAY8yNj\nzJvzPSkpP2eeH+bps/3FnoYEECT6WeCUtfY9wCeA08YY7QsQILOFP3n6AvZSipdiE5w8fYH+RLLY\n05INVOQ62IUxphaotNbOZ0//DPigtfaVda6io2c4JjE6xeOnzgHw5BMPc7C9pcgzclagb84JsiPv\nMeB1wCeNMXuBEHB1oyuMj08HGTsvotGQs+MXa+yzPxnmkbd30NRUx9mfXuTE8c6CzwHc/n9/Y/wg\ngkT/deApY8x5IA08Zq1N38XcZJvZF23mWF8b0WiI587Hij0dySFn9NbaJeDRAsxFytSxvrY1f5bS\npB1yIo5R9CKOUfQijlH0Io5R9CKOUfQijlH0Io5R9CKOUfQijlH0Io5R9CKOUfQijlH0Io5R9CKO\nUfQijlH0Io5R9CKOUfQijlH0Io5R9CKOUfQijlH0Io5R9CKOUfQijlH0Io5R9CKOUfQijlH0Io5R\n9CKOUfQijglyfHqMMW3AL4B3WWsH8jul8tWfSDI6uUB7a12xpyKyrpxbemNMNfBVYDb/0ylvZ54f\n5umz/cWehsiGgizvvwR8BbiS57mUrf5EkpOnL2AvpXgpNsHJ0xfoTySLPS2RNVV4nrfumcaYjwB7\nrbVfMMacA/4ywPJ+/RvcxhKjUzx+6hwATz7xMAfbW4o8I3FQRaAL5Yj+x0A6e/J+wAKPWGvHNrhN\nb3x8Ougkt1w0GqIY4z97Pg5AU1Mds7OLnDjeWfA5FOu+a/ySGT9Q9BvuyLPWPnTj51Vb+o2Cd9a+\naDPH+tqIRkM8dz5W7OmIrOtOXrJzctke1LG+tjV/Fik1gV6yA7DWvjOfExGRwtCbc0Qco+hFHKPo\nRRyj6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hF\nHKPoRRyj6EUco+hFHKPoRRwT+NtwReS10l6ayzNXGEzGaZ/ZwX3NR4o9pZwUvcgdSHtpXpm5ymAy\nxkAqzlBqmLnlOQB2NkT4h7cqepGylvbSXJkZZSAVYzAZZygVZzYbOcCu+h3cHz1Cb6Sbd/S8kYXi\nHdUqMEUvskraS3P1+qsMJGMMJmMM3hb5zvodvD56H73hbnoiXeyoj/jntdSHGJ8u/eoVvTjtRuSD\nyTgDqRhDyTjXl2f983fUR3j9rvvoiXTRE+5mZ0Nkg1srD4penOJ5XmZLvmq5PrN03T8/UhfmyK57\n6Il00xvuYmfDjiLONj8UvWxrnucxOjuW2fGWXa6vjjxc18qD7W+kJ9xNb6SbnfURKioCHfG5bCl6\n2VY8z+PV2TEGknEGs1vz6aUZ//xwXSvHdr+R3uxyfVfDjm0f+e1yRm+MqQS+BhggDXzcWvtyvicm\nEkQm8nEGUzFGBkd48VXL9OLNyFtrWzi2+w3+c/Jow07nIr9dkC39+wDPWnvcGPMQ8AXg/fmdlsja\nPM9jbO7aLXvXpxZv7jFvqQ3xpt330xPuojfSTbRhl/OR3y5n9NbaM8aYH2RPdgDJvM5IZBXP8xif\nu+bvXR9MxphcFXmotpkH2o4SH6yhOb2Hv/vgcUWeQ6Dn9NbatDHmm2S28B/K64zEaZnIJ/zn44Op\nOKmFSf/8UE0m8hvL9eR4Nd9/4SJXLqWAJf7x6V9y4ngnfYfK/6W1fKnwPC/whY0xbcDPgXustXPr\nXCz4DYrzPM9j7Po1fjM2wG/GBnh5bJCJuZuLyZa6Zu5t6+VIWy/3tvWyL9T+mi15YnSKx0+dA+DJ\nJx7mYHtLQe9DCQm0xMkZvTHmUWC/tfaLxpgW4Jdkol9c5yre+Hjx3pUUjYZwdfxyuO+e5zExn8y+\ndz2zNU8upPzzm2ua6Al3ZV4nj3TT3tiWc7n+7Pk4AE1NdczOLnLieOfd35lNKIHHP1D0QZb3zwBP\nGWN+nL38pzYIXuQ1JuZ+y0Aq7r9WvjrypppG7o8eyb4Zppv2pjYqK+7sE9/7os0c62sjGg3x3PnY\nVk9/2wmyI28W+OMCzEW2id/OJ7N71zOvlU/M31yuN1U3cjR6xH/v+p6m3Xcc+e2O9bWt+bOsTW/O\nkbuWnE8xkIwxMjzCi1ctE/O/9c9rrG7g6K77/OX6VkQud0fRyx1LzqcYXLVcv7Yq8obqhls+oLKv\nuV2RlxhFLzmlFiZvWa6Pz0345zVU1/O6XffQG+7mwa7X07jUoshLnKKX10gtTDLkvxkmztjcNf+8\n+qp6juy8h55I5h1v+5v3+pFHI8Xdey3BKHphcmHq5nI9FWNsdnXkdRzZ2UdPpJuecBcHQvu0JS9z\nit5BkwvTDKVi2ZfR4rw6O+afV1dVy707Db3Zj5rub95LVWVVEWcrW03RO2B6ccb/LPlgMsboqshr\nq2q5d4fxl+sHmvcp8m1O0W9D04szq5brcUavv+qfV1tVyz07ev3XyQ+G9ityxyj6bWBm8TqDqTiX\nRy7x66v9XLk+6p9XW1lDX6SH3kg3PZFuDily5yn6MjSzdD27dz2zNV8deU028hvL9YOh/VRX6n+z\n3KTfhjJwfWmWoVTcf17+ysxV/7yaympM5DA94W4e7DxCa3qnIpcN6bejBM0uzTKYGmYwlXnH25WZ\nUbzsJ5arK6v9Pes9kW4OtRygJht5sT/lJeVB0ZeA2aU5hlJxf+fb5Zmrt0R+ONzpfwqto+UANVU1\nRZ6xlDNFXwRzy3MMpYb95frl6Ss3I6+oykSe/Y63jpaDily2lKIvgLnlOWKpi/53vF1aFXlVRRVd\nrR30RrrpjXTR0XKIWkUueaTo82B+eZ6h1HB2uR5nZPrybZEf8pfrna0Hqa2qLfKMxSWKfgvML88T\nm0wwmIwx/KuLxJMjpL00AJUVlXS2HqI3+xVQXa2HFLkUlaLfhPnlBeKTFxnMvow2Mn3Zj7yqopKO\nlgP0ZN/x1tXaQZ0ilxKi6ANYWFkk7j8nj5OYvnTLlvxQaL+/XH+w+z6mU0tFnrHI+hT9GhZXFolP\nJrJ712NcnLo18oOh/f43tna3HqK+ut6/bn1NPdMoeildip6bkd9YriemLrHirQBQQQUHW/b7H1Dp\nbu24JXKRcuNk9IsrSwxPJrLveIuTmBpheVXkB0L7Mu94C3fRHe6kQZHLNuJE9EsrSwxPjdxcrk/e\nHvle//jk3eEOGqobijxjkfzZltEvrSxx0Y88zvDUCMvpZSAT+f7QXv8db92tnTTWKHJxx7aIfim9\nzMXJkcxW/KUE9lr8lsj3Ne/xl+uHw5001jQWecYixVOW0S+ll0lMXfK/GWZ48iJL2ciBTOThzKfQ\nDoc7aVLkIr6yiH45vUxi6rL/nDw+mWApffNlsX3Ne/yX0N7S/Trmp3TgXJH1lGT0y+llRqYvM5DM\nfNQ0Nnnxlsj3NrVn3wzTxeFwF821Tf55obpm5tFnykXWUxLRr6RXSExfZjC74y2WGmZxVeR7mnb7\ne9cPhzsJ1TYXcbYi5W3D6I0x1cA3gA6gFvi8tfYHdzvoSnqFkelX/G+GiU1eZHHl5tGv25t2+x9Q\n6Ql3KXKRLZRrS/8ocM1a+2FjzA7gl8AdR7+SXuHSzCsMJuPZyIdZWBX57sY2f+96T6SLltrQnQ4B\nQH8iyejkAu2tdZu6vmyeHvvykSv67wDfzf5cAcHeVL6SXuHyzBX/dfJYapj5lQX//N2NUf918sPh\nblrrNhf57c48P0xNbRV/84dHt+T2JDg99uVjw+ittbMAxpgQmfg/k+sGv/i/T/Ly2BDzK/P+f2tr\n3MUD4fv9JXtrXctdTvtW/YkkZ54fxl5KAXDy9AVOHO+k71BkS8eR19JjX34qPG/jl7eMMQeAZ4Av\nW2u/lesG/+jbn/D2NLdxb1sv97X1cG+0lx2N4S2a7voSo1M8fuocAE8+8TAH27f2D4usT499yagI\ndKGNojfG7AbOAZ+01p4LcoMTs0kvfb3wLwo8ez4OQFNTHbOzi5w43lnwOUBxv4a6WGPrsS+Z8QNF\nn6vOTwNh4LPGmM8BHvBea+3CelfY2Rhh/Hrh7/i+aDPH+tqIRkM8dz5W8PFdpse+vORc3m+C5/hf\nW+e29Bq/ZMYPtKWvzPdERKS0KHoRxyh6EccoehHHKHoRxyh6EccoehHHKHoRxyh6EccoehHHKHoR\nxyh6EccoehHHKHoRxyh6EccoehHHKHoRxyh6EccoehHHKHoRxyh6EccoehHHKHoRxyh6EccoehHH\nKHoRxyh6EccoehHHKHoRxwSK3hjzZmNMoOPTi0hpy3V8eowxTwB/Bszkfzoikm9BtvRDwAfyPRER\nKYyc0VtrvwcsF2AuIlIAOZf3mxGNhvJxsxq/xMfW+MUfP4g7ib4i6AXHx6c3MZWtEY2GnB3f5fuu\n8YP/wbmTl+y8zU1FREpJoC29tTYBvC3PcxGRAtCbc0Qco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUc\no+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj\n6EUco+hFHKPoRRyj6EUco+hFHKPoRRyj6EUco+hFHJPzqLXGmArg34CjwDzw59baeL4nJiL5EWRL\n/36gzlr7NuDTwD/nd0qyGf2JJC8OXSv2NKQMBIn+OPDfANbanwFvyuuMZFPOPD/M02f7iz0NKQNB\nom8BJledXjbGaF9AiehPJDl5+gL2UoqXYhOcPH2B/kSy2NOSEpbzOT0wBYRWna601qY3ukI0Gtro\n7LxzafxoNMTB/WEeP3UOgL/+kzdwsL2lYOOvNZ9icn38IIJE/wLwB8B/GGPeAryY6wrj49N3O69N\ni0ZDzo1/9ifDPPL2Dpqa6jj704ucON5Z0PFvcPGxL7XxgwgS/feAdxtjXsie/uhmJyX5sS/azLG+\nNqLREM+djxV7OlLickZvrfWATxRgLrJJx/ra1vxZZC3aISfiGEUv4hhFL+IYRS/iGEUv4hhFL+IY\nRS/iGEUv4hhFL+IYRS/iGEUv4hhFL+IYRS/iGEUv4hhFL+KYCs/zij0HESkgbelFHKPoRRyj6EUc\no+hFHKPoRRyj6EUcE+R77zfFGPMB4EPW2j/N1xirxiqJI+saY94MfNFa+3CBx60GvgF0ALXA5621\nPyjg+JXA1wADpIGPW2tfLtT4q+bRBvwCeJe1dqDAY18AUtmTw9bajxV4/L8HHgFqgH+z1j613mXz\nsqU3xvwL8HmgIh+3v4aiH1nXGPMEmV/8ukKPDTwKXLPW/g7w+8CXCzz++wDPWnsc+CzwhQKPf+MP\n31eB2SKMXUfm/r8z+0+hg38IeGv29/93gQMbXT5fy/sXKOwBMkrhyLpDwAeKMC7Ad8jEBpk/tEuF\nHNxaewb4i+zJDqAYR9D8EvAV4EoRxj4KNBljfmiM+VF2xVdI7wFeMsY8C3wf+M+NLnxX0RtjHjPG\nvGiM+fWqfz9grf3u3dzuJhT9yLrW2u8By4Ucc9XYs9ba68aYEPBd4DNFmEPaGPNN4F+B04Uc2xjz\nEWDMWvs/FG51udoscMpa+x4yG7vTBf792wU8AHwoO/7TG134rp7TW2u/Qea5ZLHd8ZF1txtjzAHg\nGeDL1tpvF2MO1tqPZJ9X/9wYc4+1dq5AQ38USBtj3g3cD/y7MeYRa+1YgcYfILPSw1o7aIyZAPYA\nrxRo/Ang/6y1y8CAMWbeGLPLWnttrQtvl733L5B5LkvQI+vmUcG3NMaY3cAPgb+11n6rCOM/mt2R\nBJkdqSvZfwrCWvuQtfbh7A7UXwEfLmDwAI8B/wRgjNlLZgN0tYDjPw/83qrxG8n8IVhT3vbeF1gp\nHVm3GJ9g+jQQBj5rjPlcdg7vtdYuFGj8Z4CnjDE/JvM79Slr7WKBxr5dMR7/r5O5/+fJvHrxWCFX\nmtba/zLGvMMY83MyG52/yh54dk36lJ2IY7bL8l5EAlL0Io5R9CKOUfQijlH0Io5R9CKOUfQijlH0\nIo75f3LbxP56ayccAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108fa0290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-1, 5, 1000)\n",
    "y = coefs[0]*x + coefs[1]\n",
    "\n",
    "plt.plot(A[:, 0], b, '*')\n",
    "plt.plot(x, y)\n",
    "plt.xlim(-1., 6)\n",
    "plt.ylim(-0.5, 5.5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are not sure about the result. Just check it with another method. For instance, I double-checked with R:\n",
    "\n",
    "```r\n",
    "a <- data.frame(x=c(0, 1, 2, 3, 3, 4),\n",
    "                y=c(2, 4, 0, 2, 5, 3))\n",
    "\n",
    "ggplot(data=a, aes(x=x, y=y)) +\n",
    "  geom_point() +\n",
    "  stat_smooth(method = \"lm\", col = \"red\") +\n",
    "  xlim(-1, 5) +\n",
    "  ylim(-1, 6)\n",
    "```\n",
    "\n",
    "outputs:\n",
    "\n",
    "<img src=\"images/linear-regression-r.png\" width=\"300\" alt=\"Fitting a line with another method (in R)\" title=\"Fitting with R\">\n",
    "<em>Just checking with another method</em>\n",
    "\n",
    "You can also do the fit with the Numpy `polyfit()` to check the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.21538462],\n",
       "       [ 2.2       ]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.polyfit(A[:, 0], b, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's good! We have seen how to use the pseudoinverse in order to solve a simple regression problem. Let's see with a more realistic case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 4.\n",
    "\n",
    "To see the process with more data points we can generate data (see [this nice blog post](https://mec560sbu.github.io/2016/08/29/Least_SQ_Fitting/) for other methods of fitting).\n",
    "\n",
    "We will generate a column vector (see `reshape()` bellow) containing 100 points with random $x$ values and pseudo-random $y$ values. The function `seed()` from the [Numpy.random package](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html) is used to freeze the randomisation and be able to reproduce the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "x = 5*np.random.rand(100)\n",
    "y = 2*x + 1 + np.random.randn(100)\n",
    "\n",
    "x = x.reshape(100, 1)\n",
    "y = y.reshape(100, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We will create the matrix $\\bs{A}$ from $\\bs{x}$ by adding a column of ones exactly like we did in the example 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.48234593,  1.        ],\n",
       "       [ 1.43069667,  1.        ],\n",
       "       [ 1.13425727,  1.        ],\n",
       "       [ 2.75657385,  1.        ],\n",
       "       [ 3.59734485,  1.        ],\n",
       "       [ 2.1155323 ,  1.        ],\n",
       "       [ 4.90382099,  1.        ],\n",
       "       [ 3.42414869,  1.        ],\n",
       "       [ 2.40465951,  1.        ],\n",
       "       [ 1.96058759,  1.        ]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.hstack((x, np.ones(np.shape(x))))\n",
    "A[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now find the pseudoinverse of $\\bs{A}$ and calculate the coefficients of the regression line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.9461907 ],\n",
       "       [ 1.16994745]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_plus = np.linalg.pinv(A)\n",
    "coefs = A_plus.dot(y)\n",
    "coefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We can finally draw the point and the regression line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAECCAYAAAAVT9lQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0lGWe6PFvVfY9IQshgSSsTwggyKYCyqKCIm6tGEAa\n7GWcsZ3b0zM9ntveuT09c+b0TDtLt9PTbc89PdqiDc0iCm4sKqiAtCAgIkketpCENQtZyVrL/SNJ\ndRKqUpXKW1v4fc7xdFWo931/qc77e5/9MdntdoQQNzdzoAMQQgSeJAIhhCQCIYQkAiEEkgiEEEgi\nEEIA4Z5+UCl1G/AzrfVCpdQ04JeABWgD1mitq3wUoxDCxzwqESilngN+C0R1/ehF4Fmt9SLgLeBH\nvglPCOEPnlYNzgCP9nhfqLU+0fU6HGgxNCohhF95lAi01m/RWQ3ofn8VQCk1B3gW+IVPohNC+IXX\njYVKqULgJWCp1rrGuJCEEP7mcWNhT0qp1cDTwAKtdZ0nx9jtdrvJZPLmckKIwXF745k8nXSklMoF\n/gDMA6qAMqAesAOfaK3/0c0p7FVVjR5dKxikpycQSvGCxOwPoRYvQHp6gttE4HGJQGtdBszpepvq\nbVBCiOAjA4qEEJIIhBCSCIQQSCIQQiCJQAiBJAIhQkZJWS0lZbU+ObdXA4qEEP63fX8pAPm5KYaf\nW0oEQgS5krJaXlh/FF1Rh66o44X1R/stGXhTcpASgRBBLj83hYTYCH788iEAVi9RZKfFufx8d8mh\nW3p6gttrSCIQIgQcLqnkobl5AHxRUkn2vNE3fKakrJbt+0vRFZ3Tf/7rza9IS4rhzpk5bs8vVQNx\n0/Flo5uvZKfH88idY3jkzjFkuSgN5OemsHrxBMf7ljYrFZVNPPjD7R+7O7+UCMRNx5eNbr4yKz/D\n6eu+uksODc3tfHzsUvePv+fu/FIiEDeNgTa6haLukkNibCQzVXp3dWK5u+MkEYibRt+i8+olKqRK\nBZ7oLi1kp8ezZtkY2jKOAxS5O06qBuKm4kmj21AQm17LTz//fzS0N/LOf/xmi7vPSyIQN5Xs9HjH\nU/NwSWWAozFem7WdbWfe49OLBwkzhfHw2Ps9Ok4SgbipeNroFopK68t5rWgjlS3VZMVlsrZgBSMT\nsjw6VhKBECHOarOy8/xH7Czbg91u5+5Rd/HgmCVEhEV4fA5JBEKEsKvXK1lXtImyxgpSopJZkLqU\nkZG5A0oCIIlAiJBkt9vZd/Egb555jw5bB7Mzp/PEhIf55aZijlA64N4QSQRChJi6tnrWF79B0TVN\nXHgsawoKiW0ZxS83FTuGF7+w/igPzxvtcUKQRCBECDla+RUbS97kuqWZicMmsHricpKjkgAGNDGp\nL0kEQvhR90jGgRbdWywtbD61nUNXjhJhjqBwwiPcmX0HPTcNGswYCUkEQviRN/McTtWe5bWiTdS2\n1ZGbMIq1BYUMj7ux63MwYyQ8TgRKqduAn2mtFyqlxgKvAjbga631swO6qhBByNuntafn7jlF2JM6\nfIe1g3fO7WJPxT5MJhNL8+7hvry7CTOHOf38YMZIeDTXQCn1HPBbIKrrRz8H/o/Wej5gVko9PKCr\nCuEjg5livH1/6Q2LehhloPMcLjRe4l+/+C8+qviU9JhUfjjjezwwZrHLJDBYnpYIzgCPAq93vZ+h\ntd7X9XoHcC+w3eDYhBgwb4re3jytveFJHd5mt/Fh+Se8e243VruVu7Lv4JFxDxAVFmloLH15lAi0\n1m91bYLareemio1AkqFRCTFAg7mZB7oUmLfc1eGrW67xWtFGztafJzEygdUTlzMpNd/wOJzxtrHQ\n1uN1AuDR1uierJ0WTEItXrh5Y05PTyBnZDJ/+W97Afj+ilvJyUz0+PgPjl5k5WIFQHF5HdMmZvZ7\nLW8s7XFcz9d2u52PSw/yu2ObabW0cdvIW/mzmatIjIr36jre8DYRHFVK3aW1/hS4H9jjyUGhtJ10\niG5/fVPHvPuzUkfRe/fB8zw8gO6z5NiIXk9rVzEZ/R03tjfxh5KtHK8+SXRYNGsmFjI7czptDXaq\nMOY6vly89G+B3yqlIoBi4A0vzyOEYQbTfebvWYklZbWUXj/NJzU7aexoYnzyGL45sZDUmMAslGKy\n2+3+upY9lJ5WN/vT1V9CLWYj4m21tPGTHa/SFHcWM2HckTqfFbcsxmzyzYJh6ekJJnefkQFFQvjR\nXn2CbWVvYYlrwtacABXT0JEpmKcGdtVAWbNQCD+w2qy8c24XWy/+Hkt4Ex2XRtN28g5a6uOoqLru\ndCFVfy67LiUCIQzWd4TiletXWVe0kfLGiwyLTmFMx52k5GbRmN7O3q4lx511Wfpz2XVJBEIYrPsG\nnpCTxKcXDrLt7Ht02CzcnjmTxyc8xNdnGpiVn8G2feeYqdLJSovrNcDIXwOcepJEIIRBet3AEa08\nt/MXtEZdJT4ijqcmrWJa+mQAZuVHA3/agwB693L4a4BTT5IIhDBI9w38D29tJyLvJK3hFian5rMq\nfzlJUTf25ffXZenvZdclEQhhkOaOZl75egOR485gJpz8sPn8xS1Le60Z4Cl/L7suiUAIA5RcO83r\nxZups9czOjGHNQUrKCu3eZUEwP8DnCQRCDEI7dYO3j67g70X9mM2mVk2egmLcxcQZg4jwz/zhQwh\niUAIL5U3XmDdyY1caa5keGwGawsKyU0cFeiwvCKJQAScL1cG8gWrzcoH5Z/wXulubHYb80fO5ZGx\n9xPp4zUDfEkSgQg4IwbO+CuZXGmq4sVjL3OuvoykyES+OfEJTE3pnLtwnfxc14kg2JOdJAIRMCVl\ntfx8y3FDBs74ehSe3W7ns8uH2HrmXdosbczImEqhepS4iFhe2HnU7bX9OUrQGzL70IVQmxUHoRlz\ns9XuWEzkn75724AHzvQdhadGJRs+Cq+hvZENJW9worqY2IgYCsc/wszMWz26tj/ic0dmH4qgd+D4\npUENnPH1KLzjVSfZUPIGTR3XmZAyjh/M+xb26xEeXzsQowS9ISUCF0Lx6RqKMetLDaisziXFDpdU\netVnvm3fOcdrk8k0oJWJXGm1tLL19Dt8dvkw4eZwHh57PwtGzmV4RlKv79iTa/sivoGQEoEIevOm\nZjtuLG8Hzhg9Cu9s3XnWFW2kpvUao+KzWFOwgqx452sYenJtf48S9IaUCFwIxaerxDw4FpuF90o/\n4IOyjwFYnLuQpaPvIdz8p+dlMMXrKSkRCOGhS01XWFe0kQtNl0iLHsaaghWMTc4LdFh+I4lADHn9\n9eHb7DY+rtjP9nM7sdgszBkxm8fGLyM6PNrfYQaUJALhd31vTF8PtnHVh3+ttZbXi7dwqvYM8RFx\nPDnpSW5Jn+STGIKdJALhd31vTF8NtnG10o/KSebw1WNsPrWNFksrU9IKeDL/cRIi/behSLCRRCD8\npu+N+ZNXDhEWbub8pQbA+CW5nPXhJyeZeOXkeo5WfkVUWCRP5j/OHSNmeT1deKiQRCD8pu+N+fRD\nk0hJiXWMLPTFYJueK/3sLDrCWfOn1Lc3MCYpj7UFhaTFpBp6vVDlVSJQSoUD64A8wAL8mdb6lIFx\niSGq7xJcsbGRPl2SKzs9nqnjk9l29n2OXvgMs8nMQ2Pu497cBT7bUCQUeVsiWAqEaa3nKqXuAf4Z\neNy4sMRQ1XdwTWJiNCor2/F+IDxpZMzIauNnh/+Tq81VDItI4+lpTzIqIdvL6I2LK9h4mwhOAeFK\nKROdW6K3GxeSGMr6LsHVc4DOQEcW9tfIaLVZ2V22l/fPf4jNbiOhaQKxjbcw6k7fJgF3cQUrbxNB\nEzAaKAFSgWWGRSSEG+7W/a9srmJd0SbON5QTH5ZAbOUMys5GU0kTP3nlECvvHu+Tm9RZY6ivrmU0\nbxPBXwM7tdZ/p5TKBvYqpSZrrfstGXi7r3yghFq8cHPEnJ6eQM7IZEcj4/dX3EpOZiJ2u50Pz+7n\ntS/foM3azrzc2XxneiHV1yyOz0ZFhXHnzByfxNs3LiOu5S/eJoJrQEfX67qu84S5OyiUxmiH6Jjy\nIROzu3r27s9KHY2Muw+eZ8GsNNaXbOFkTQmx4TF8e9IqZgyfRnO9lQ07igkPM2Gx2jlTUc/fvviJ\n192U7r7jDTuKSU2MoqahbdDXMoonidbbRPAi8IpS6lMgAnhea93i5bmEcOhOAO7q2T0bHd84tp+f\nHlrH9Y5m8lPG882CJ0iOSnKc7/zlBizWP02umztlhM9uzKnj0rj/tpygX3+gL5l96MJQeroaxajW\n8J7n6RvzT145RHV9Cy1tVqD/FX1aLC1sOfU2n185QoQ5nDuGLeSWxBlMzBvW6xoXq5ocN+bCW7NI\njItyuiaAJ7+fJ99xoNcf6EtmHwpDGdUa7uw83Q1tFZVNvT7r6ol6uvYcrxVv4lprLRlRmdyb8RC7\nPq2hiNP843du63WNwyWVjs1GTSYTWS6e0Eb9fqGw/kBfUiJwQUoEf2LUunvOzrN22SQyk6IAbnhy\nJ8RG3vBE7bBZePfcLj4q/xSA+/IWcWhvIjUN7Y5SRExUWK8ShcpJ7rXZaN9uyoH8fiH6dyElAjF4\nRq275+w8U8alOW6swyWVpCZG0dZhcxTfez5RLzZd5tWTf+DS9SskRaRwS/jdnDwYzoWqul7XSUmI\npqXteq9Yu4v9zsYqhMq6gr4kiUB4xKjdefueZ9rEziXASspq+aKkkpqGts7PFV9FjUpmVn4GNruN\nPRX7eOfsTix2K/OybuP80ZGU2cNZvXiC4wbuFhluviFWd8V+f+8+HGwkEQiPGFXvdXWe/NwUnnlk\nsuOmfubRKWSnxVHTUsvrxZs4XXeOhMh45g+7n6+OhHG6qxj/4pbjzJk8nNLLjSTGRqJykjlZes1R\nFdi27xwvrD/qdu+EUKzXG0naCFwI0bpgSMfc3dpeXd/ZEz1lRiubT22n1drK1PTJrFTfICEyvldb\nQuGicSyZneO4eWflZ9zQDtDz897sneAq3lAhbQQiILztZux+Kv/zhoNUxn/OseJLRIdFsXriE9ye\nOcOxZkDPYnxre2ejYH/biN/sxX5PSCIQhvO2Gy4hJoJ/eONdKod9jimyjai2dArHP87sEb1v3IEW\n42/2Yr8nJBEIw7ibDNSfNms7x1r3UjXsj2Az0VGu+NEDTzDKyfDY/p7+zgz08zcjSQTCMO664VxV\nGUrry3mtaCOVLdXEM4zJkYtIGJvGUV3tNBEI40kiEIbqrz7et8pgtVnZ/PU7vFm0E7vdzt05dzGi\nYzp3TMyipKyWkvJav8d/s5JEMMT5e7UcZ/VxZ1WGu2Yn8mnd+5Q3XiAlKpk1BYVMSBnrOE930hD+\nIYlgiPP3ajnO6uO9qwx28mfUs+niVjpsHdyVdxsP5SwlJjwGGFw7Q/fx3dcUnpNEMEQN9oYy2uGS\nShbPSafI+jG7L1UQFx7LmoJClkya26tffrDDfUNxmbBgIIlgiAq28fMd8Rc52vgR1+3NZEeN5tmZ\nT5IUlej0s970+wdb4gs1kgiGsGAYSNPc0cLmU9s5XHuUCHMEhRMe5c7s2/vdUMSbfv9gS3yhRhLB\nEObtQBpv6tnOjjlVe4bXijZT21ZHbuIo1hasYHhsuttzedvvHwyJL1RJIhjCvL2h+tazXSWGnj/v\neUyHtYO3z+1kT8U+zCYzD4y+lyW5iwgzu13WclBkBKH3JBEIB1f1bFcNcNv3l9Lc2kFsdITjmH/a\n9BHtWUe41lFNRmwaTxWsJDdxlF/ilxGE3pNEMIQNtIjft549d0qm08QA9Pp5XmYCYCd8RCmVGWew\nddi4K3sOj45bSmRYpMG/lfAFSQRDmDddaT3r2TUNbb0W/ujZANczYYwYYaIt9zh19itE2GP5ztSV\nTEpVBv4mwtckEQxBg+lK61vPdtUAd7ikkgfn5HLJrvnK8iE2u4Vb06egwu5iUmqub34x4TOSCIag\nwXSlOatnO2uAGzbMRLF1H0XVJ4kwR/LNiSuYNfzWG7oF+6uelJTVcqW+zbF4qQgcSQRDlFFdac4S\nw4nqIt6/9gaNHU2MTx7DmoJChkU7L230Vz3Zvr+UiMgw/mb5VK9iE8bxOhEopX4EPETnTkcvaa1/\nZ1hUYtB80ZXWamnjzTPvcODSIcJNYXxj3DIWjpqH2WS+4bP9VU9kFGDw8SoRKKXmA3dorecopeKA\nHxoblhgsZ0/ywUzIOVd/nnUnN1Ldeo3s+BE8VbCSrPhMl5/vr3oiowCDj7clgiXA10qpbUAC8Jxx\nIQlf8aYXwWKzsKP0Q3aVde7we2/OAh4Ys5gIs/s/nf6qJ93/FhcXJaMAg4C3iSANyAGWAWOAt4F8\no4ISxvK2KH75+lXWFW2kovEiqdEprClYwbhkz2/Y/qon3f+Wnp7A+/vODvA3EkbzajlzpdS/AJVa\n6190vf8SuEdrXd3PYX5bN13cqOxKA3/5b51P9V8/t5CcTOcz/wBsdhs7T3/M+q+20WHtYOHoOay9\n9XFiI2L8Fa4wls+WM98PfB/4hVIqC4gFatwdFErrwYfo+vUuY979WamjmL774HmXO/TWttbx++It\nlNSeJj4ijqcKVjItfTLX6yxcx/jvI9S+51CLFzpjdserRKC1fk8pdadS6hCd2eZ7Wmt54hvE29l/\nPfvk+57Dk16EL65+yUb9Fi2WFianTmRV/uMkRcnioTcDr7sPtdY/MjIQ8Seezv7re0zPPvm+5+hv\nQk5zRzMb9VscqTxOZFgkq9RjzMma7XLNAFkObOiRAUVBZKCz/5wd85OXPweTiYrKpl7ncHXTllw7\nzevFm6lrq2d0Yg5rClaQEZvWb5yyHNjQc+NIEBEw+bkprF48wfG+5+w/XVHHC+uPOp7Gro55+uHJ\nPP1ggeP96iXK6Q3bbu1gy6nt/NeXv6WhvZEHxyzhr6c/028SKCmrdWwo6iqekrLaG34mgp+UCIKM\np7P/nB3T3Sdvt9v7HV5c3nCBdUUbudJcyfDYDJ4qWEFO4ki3sXkyEEhKC6FJdkN2IVCtwz138j1c\nUsnFqibHv5lMJqet/d3H9OyT7z7Htn3nyM9JIT83BavNygflH/Ne6QfY7DYWjJzLw2OXEhkW4TQW\nZ20B3TsW942nbxVFjUr2aKxCqLXCh1q8ILshhyRPZ/+5O6abLq9Dl9eRmm5jXdFGShvKSI5KYvXE\n5UwcNqHvqXpx9nR31fsgw4ZDm5QIXAjRzO+I+U9P6FrC0i8Qlauxmy3MyJjKCvUosRGxLs/j7dN9\n275zVNe3dsaSHONyrIKrmENBqMULUiK4qeXnpmCKaOMXBz8iLKWKqLBoVqrlzMy81aNjvXm6Z6fH\no8s7k8fUcf33PIjgIolgiDpe9TW/K9lMWEorw0zZTDYvYmbmFI+PH+h6BiVltew5csFRigBIiImQ\nRsMQIYlgiGmxtLL19DscvHyYMMJ5fPxDzB85hyP9TgO50UDXM5A2gtAmiWAIKa46zS8P/Y6a1lpG\nxWexdtJKRsQNp6SsloQY5z0DrnizNLhsMBK6JBEEOU+G81psFt4r/YAPyj4G4L7cRdw/+h7Cu9YM\n8FffvmwwErokEQSIp+P13d3El5qu8GrRH7jYdJnhcWmszn+CMUl5jmv4c0kw2WAkdEkiCBB38wfK\nrzZy7HS1y5vYZrext2I/28/uwGq3MjdrNn9++0oa6zoc55F6u/CUJAI/c/aUvnV8GjnDExw3eXeS\ncDW8+FprLa8XbeZU3VnM1ijS6uewatGDREdE00hHr+tJvV14QhKBnzl7Sv9+l+bY6c5W/Z5J4sUt\nx5kzOZO0pGi+KKkka24eh68eY5PeRqu1lZjWbK4VTeC6JYoX1h9l7bJJN+wRMNB6u0wxvjlJIgiA\n7qd0dX0rL24+Tk1D52i87ftLmTsl05EI7pk5itzhnQuDVF1v4OWT6zlW+RVRYZE8mb+cnPB8/v6r\nw0BnQpkyLu2GUW8DrbfLpKGbkySCAOj5lN51qJxNe84AnTfz4eKrjqJ8a7uV7ftLaYm8TPuIo9S3\nNzImKY+1BYWkxaSybd+5XsX+aRNdLy/ujuw1cHOTRBAAPZ/MLW2W3nX4riRRUlbL6x+cpDruS8JT\ny6HNzB2pC1g19T7HhiJGdtdJw+LNTRLBAPii/tz3Zu5+HZPShDn/AOHtNdia4/mzqU8yfdTYXsca\n3V3nzbBikGrEUCCJYACMqj/vPlQOwOLZOTfczFablV1le9hx/iNsdhs55lsYlzibijIz00e5PqcR\nG4oOtIQh7QlDhyQCDxhdf95+oPMGWjw7p9fPrzZXsa5oI2UNFaREJTM7fjEPTZ0JuL8xjdhQ1NMS\nhrQnDD2yHoELfeedX6xqctSf/+m7t3lVf959qJztB0ppabMCEBMVxu0Fw5mpMqiO0Lx5+l3abR3M\nGj6dJyY87NGGIt6uHTBYRnwfEHrz+0MtXvBsPQJJBC70/T/c1RJdA3X0VCW/evNrAP7XY1PY8cUp\napIO0Rp9mdjwGFaobzBj+MCe6kbdlANh1PcRajdWqMULsjCJoYxqod99uIIJo5JobrXwP/v2YM/+\nClNEB9FtmRTmPs6M4TnuT9JHIDYUlQlGQ8ugEoFSKgP4gs59D08ZE5LvDKaV26gW+rzMRB6eP5It\np96m5soRsJppPz+R5x9czsj0eK/O2XdDUX+05ssEo6HF60SglAoH/htoNi4c3wpUK3fPG3PGdDM/\n/fwX1LbVkWhKZ3LUIuLGp3BEV3mdCPrelC+sP+q4nhCeGEyJ4N+B3wDPGxSLzwS6lXv7/lLsWJnQ\ncZU95fswmUzcn3cP6W1TuG3iCMCY4vWJM9Wse/ektOaLAfNqpyOl1FN0bov+AR5suWy0ge6m03c3\nIFe7/xite2egU9UVlCXt4KPyT0mMSOZvpj/DsjGLHUkAjCleTxmXFpDfU4Q+r3oNlFKfALaut9MA\nDTykte7vsWZY98TzL+0H4F++N8/jYzbsKnG8NgErl+QbFQ4nznTOHJzSZ+Vem83G60fe490zOzCZ\n7dwx4naembOC6HDvB/2448vfU4Qs3/QaaK3nd79WSu0F/txNEgAYdLdL3yL+3774icdF3+TYiF6t\n3O5iGUg30bp3TwLwv5+c7vhZTcs1XivexJm6UqLMMRSELSChJo/G2nYaaffovAOVnp4w4N8z0EKt\nOy7U4oXOmN0xovvQbwMRPJkY46rF3Bet3M7aHh6am0d91Dm2nNpOq7WN3OjxfG/mKuIj4/zSzSat\n+cIbg04EWutFRgTiKXcTY/zdM9Bz/YBv3D2Kj6ve5ctzXxMdFsWaiYXMzpyOydRZMpMbUwSrkBtQ\n5Gogi796BnqWOLbvL6W6voWH5uZRZSvjpaJf0U4L45JHs2ZiIakxwwy9trt4hPBWyCUCV0Vfo+bT\nd99YrupV2/eX0tzaQWx0RGfSMVs4WPchTXFnMBPGo+MeYNGoOx1rBviazAAURgi5RNAfIxbq7L6x\n7pzZe6hv3xJHXmYCprg6Isd+RVN0M1lxmTw1aSXZ8SNuOOdAefKUd1YCcrZmoRCeGFKJYDDj3/ve\nWM+/tJ+ls3McN2OvEofJRnj2aWLCjmPHTp55Gj+Y9QQRZmO+Tk+e8s5KQM7WLBTCEyE7+9AXdeOe\ns/h+/dxCYsJ6d79u23eO6/ZaTlj30GivIiUqmbUFhdRdTTCkIXCgU4r7zgD87qO3uEwEwdqWEGrd\ncaEWLwzx2Ye+qBv3rFocOH6Je6ZnO/7NbrfTEHuKo42f0mG3MDZmEs/MeoKY8BgwKISBtnMMpAQk\nbQmiPyFXIvDlQhw91wzUlxpQWYkA1LXV8/viLRRfO0VcRCwr1WPcmuH5FuP96fukHsw8f2dPq0At\nXOKpUHvChlq8MERLBL5cbbdn8X7e1Gyqqho5cvU4G/WbNFtaKEhVrM5fTlJUoiHXgxuf1EbP85fV\niYUnQq5EAMatjtOf2KQwXvrs9xy+eoxIcwTfGL+MeVm3OwYHDZYvntSunlb++L68FWpP2FCLF4Zo\niQB8vzqOvnaG9Qe3UNNSS15iDmsLCsmITTf0Gv58UstqQsKdkCwR+EqHtYO3z+1kT8U+zCYzS/Pu\nYXHuQsLMYT65ntFP6hB9WoVUzKEWL4RYiSDQ3VsVjRd5tWgjV65fZXhsOn8159sk2VJdft6IeOVJ\nLYJF0CSCQHVvFZ2vYc+FT9Dth7BhY/7IOTwydinZqan9Zn4j4pWZgiJYBDwRBHIZseqWGv6n5GXa\nIqsJs0bz7PQnKUhVQRuvEL7in5kx/QjEMmJ2u52tJz7mHw78B22R1VhqMmk6PodN2+rcLoEWqGXP\nhPClgJcIwJjJQp5qbG9iQ8lWvqo+SVR4FI2nC7DWjABMPP3wZI9a7v0ZrxD+EBSJwF+NZl9VnWRD\nyVYaO5qYkDyWjKbbKYu1QCyYMHl8U3sTb6AbQ4XoT1AkAl83mrVaWtl6+l0+u3yIcHM4j41bxoJR\n8ziiqxmX0XnNgSQgb+KVsf4imAVFIvClc/XnWXdyI9Wt18iOH8FTBSvJis8E/NNq78+Vkwa7Lbq4\neQ3ZRGCxWXi/9EN2l+0FYHHuQh4YfS/hBq0Z4Cl/jSA0Ylt0cfMakong8vWrrDv5ByqaLpEancKa\nghWMSw5cg54vGxelO1MYYUglApvdxicXPmPb2fex2CzMGTGLx8Y/SHR4dEDj8mVjqMwuFEYYMomg\ntrWO14s3o2vPEB8Rx6pJTzI1fVKgwwJ83xYRiG3RxdAyJBLB4SvH2HTqLVosrUxOnciTEx8nMdL9\n7i5DRd9t0YUYKK8SQdeW6K8AeUAk8FOt9TsGxuWR6x3NbNJvcaTyOJFhkazKf4w5I2YbtmZAqJA5\nC2KwvC0RrAaqtdZrlFLDgGOAXxNB8bVTvF60mfr2BkYn5rK2YAXpsa5nCwohXPM2EWwGtnS9NgEd\nxoTjXru1nW1nd/DJhQOYTWYeHHMf9+bMd7tmgIzsE8I1b3dDbgZQSiXQmRD+zsigXClvuMCrRRu5\n2lxJZmwGayetICdhpEfHysg+IVzzeoUipdQo4E3gV1rrdR4c4vVSSFablW3Fu3jj5HtY7TaWjl/I\nqlseITJoYDH7AAAHIklEQVQ80u2xJ85Us2F3CV+frQFg8thUVi3OZ8q4NG/DESLU+GaFIqXUcGAX\n8KzWeq+nx3mzxFNlczWvFW2ktKGc5KgkvjnxCfKHjae+tg1oc3t8ZlIUhQvGOhJB4cJxZCZFuY0l\nRJekkph9LNTiBdf7ePbkbRvB80Ay8GOl1N/T+bS/X2vt/s70kN1u58Clz9l6+h3abR3MHD6NwgmP\nEBsRO+Bz9RzZt/OPZcydMkKqCEL04G0bwQ+AHxgci0N9WyMbSrbwdU0JMeExfGvicmYOn+b1+XqO\n7Pu/v/0j1fWtkgiE6CHoBhR9WfU1fyjZSlPHdfJTxrN64nJSopM9OtZVz8Cs/AzHmPxLNc1Q0yxj\n8oXoIWgSQYullTdOv80fL39BhDmc5eMf5q6Rd2A2eb6aWn89AzImXwjXgiIRnKkr5bWijdS01jIq\nIZunClaQGTfc4+M9nYEnS4wJ4VxAE0GHzcJ753bzYfknANyXdzf359094DUDPH3ayz4CQjjn90TQ\nXY9PSG1lXdFGLjZdJi0mlbUFKxiTlOv1eT152suYfCGc83si2Lb/HI1xmqakE1jsVuZm3cY3xi0j\nOnxwS2zJ014I7/ktEZw4U83LOw5REXeAsMRrmC1RLMt+lPvzZxtyfnnaC+E9vyWC2vCz1GTuJszW\nhvVaBn81bw0TMuWGFSIY+C0R/OrzVwkjgoKwBWRlKIrPXGdCpr+uLoToj98SQX7aWCaZF7Jocud2\nYVKPFyJ4+C0R/OOiH1Jd3eR4L/V4IYKH3zZBvdmWDxMilAR8N2QhROBJIhBCSCIQQkgiEEIgiUAI\ngSQCIQSSCIQQSCIQQiCJQAiBJAIhBJIIhBBIIhBCMIjZh0opE/ASMBVoBb6rtT5nVGBCCP8ZTIng\nESBKaz2Hzi3Qfm5MSEIIfxtMIpgH7ATQWn8OzDQkon6UlNU6VkEWQhhnMIkgEajv8d6ilPJpm8P2\n/aWO3YyEEMYZzApFDUDP/ZbNWmtbfwd4sj2zMyfOVLNhd4ljJ6OfbznOqsX5TBmX5tX5POVtvIEk\nMfteqMXricEkggPAMuANpdTtwAl3B3i7r3xmUhSFC8by9dkaAAoXjiMzKcqn+9Snpyf49Py+IDH7\nXqjFC54lrsEkgreAe5VSB7ref2sQ53JL9i0Uwne8TgRaazvwjIGx9NJ3i3PZyUgI3wmK3ZCd6bvF\nuexkJITvBN3IwpKyWl5YfxRdUYeuqOOF9Uely1AIHwu6RJCfm8LqxRMc71cvUY5SgRDCN4KyaiAN\ng0L4V1AmAmkYFMK/gq5qANIwKIS/BWUiEEL4lyQCIYQkAiGEJAIhBJIIhBBIIhBCIIlACIEkAiEE\nkgiEEEgiEEIgiUAIgSQCIQSSCIQQSCIQQiCJQAiBJAIhBJIIhBB4uVSZUioR+D2d+x9GAD/UWv/R\nyMCEEP7jbYngb4APtdYL6Nzh6NeGRSSE8DtvFy/9OdDW9ToCaDEmHCFEILhNBEqpbwN/DdgBU9f/\nfktrfUQplQm8Dnzfp1EKIXzKbSLQWr8CvNL350qpKcAGOtsH9vsgNiGEn5jsdvuAD1JKFQBbgSe0\n1m63QxdCBDdvE8E24BbgPJ3VhTqt9aPGhiaE8BevEoEQYmiRAUVCCEkEQghJBEIIJBEIIfDxtuhK\nKRPwEjAVaAW+q7U+58trGkUpdRvwM631wkDH0h+lVDid4zzygEjgp1rrdwIalBtKKTPwW0ABNuAv\ntNZFgY3KM0qpDOAL4B6t9alAx+OOUuooUNf1tlRr/R1nn/NpIgAeAaK01nO6bqyfd/0sqCmlngO+\nCTQFOhYPrAaqtdZrlFLDgGNAUCcC4EHArrWep5SaD/wzofF3EQ78N9Ac6Fg8oZSKovN7XuTus76u\nGswDdgJorT8HZvr4ekY5A4TKuIjNwI+7XpuAjgDG4hGt9Xbg6a63eUBt4KIZkH8HfgNcCnQgHpoK\nxCmldimlPux6GDvl60SQCNT3eG/pKhYGNa31W4Al0HF4QmvdrLW+rpRKALYAfxfomDyhtbYppV4F\n/hNYH+Bw3FJKPQVUaq0/oDPhhoJm4N+01kuAZ4D1ru4/X9+UDUBCz+tprW0+vuZNRyk1CtgDrNNa\nbwp0PJ7SWj8FTAD+RykVE+Bw3PkWcK9Sai8wDXitq70gmJ2iK8lqrU8DNcAIZx/0dRvBAWAZ8IZS\n6nYg1OYlBH3mV0oNB3YBz2qt9wY6Hk8opVYDI7XWP6OzEdna9V/Q0lrP737dlQz+XGtdGcCQPPFt\nYArwrFIqi86H8mVnH/R1IniLzix6oOv9t3x8PaOFwvjr54Fk4MdKqb+nM+b7tdZt/R8WUG8Cv1NK\nfULn3+Bfaa3bAxzTQITC3wXAy3R+z/vo7J35tqsSucw1EELIgCIhhCQCIQSSCIQQSCIQQiCJQAiB\nJAIhBJIIhBBIIhBCAP8fnA2cE0xPR+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a3c3b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_line = np.linspace(0, 5, 1000)\n",
    "y_line = coefs[0]*x_line + coefs[1]\n",
    "\n",
    "plt.plot(x, y, '*')\n",
    "plt.plot(x_line, y_line)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good!\n",
    "\n",
    "# Conclusion\n",
    "\n",
    "You can see that the pseudoinverse can be very useful for this kind of problems! The series is not completely finished since we still have 3 chapters to cover. However, we have done the hardest part! We will now see two very light chapters before going to a nice example using all the linear algebra we have learn: the PCA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span class='notes'>\n",
    "    Feel free to drop me an email or a comment. The syllabus of this series can be found [in the introduction post](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-Introduction/). All the notebooks can be found on [Github](https://github.com/hadrienj/deepLearningBook-Notes).\n",
    "</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# References\n",
    "\n",
    "## Intuition\n",
    "\n",
    "- [Sean Owen - Pseudoinverse intuition](https://www.quora.com/What-is-the-intuition-behind-pseudo-inverse-of-a-matrix)\n",
    "\n",
    "## Numpy\n",
    "\n",
    "- [Numpy - linalg.pinv](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.pinv.html)\n",
    "\n",
    "- [Numpy random seed](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html)\n",
    "\n",
    "## Systems of equations\n",
    "\n",
    "- [Using the Moore-Penrose Pseudoinverse to Solve Linear Equations](https://www.youtube.com/watch?v=5bxsxM2UTb4)\n",
    "\n",
    "- [Overdetermined systems](https://en.wikipedia.org/wiki/Overdetermined_system)\n",
    "\n",
    "- [Least square solution in the triangle center](https://math.stackexchange.com/questions/471812/is-the-least-squares-solution-to-an-overdetermined-system-a-triangle-center)\n",
    "\n",
    "- [Symmedian](https://en.wikipedia.org/wiki/Symmedian)\n",
    "\n",
    "## Least square fit\n",
    "\n",
    "- [Least square fitting](https://mec560sbu.github.io/2016/08/29/Least_SQ_Fitting/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
